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Abstract 

This  research  effort  addresses  modeling  of  the  transportation  of  air  pollution  in 
the  atmosphere  and  the  numerical  analysis  of  the  partial  differential  equations  u<5ed 
in  such  modeling.  Three  Gaussian  models  are  examined  and  compared  using  example 
problems.  Several  finite  difference  schemes  are  developed  to  solve  the  partial  differ¬ 
ential  equations  used  in  air  pollution  transport  modeling.  This  study  examines  three 
Gaussian  models:  SCREEN,  AFTOX,  and  the  program  GAUSPLUM.  The  model 
GAUSPLUM  is  developed  in  this  study  and  uses  the  Ada  programming  language 
and  the  analytic  solution  to  the  advection-diffusion  equation.  Numerical  analysis 
of  several  of  the  partial  differential  equations  (PDE)  used  in  air  pollution  modeling 
is  also  examined.  The  equations  are  generally  parabolic  or  hyperbolic  PDE’s.  The 
following  are  examined  in  this  research:  the  advection  equation;  the  one-,  two-,  and 
three-dimensional  advection-diffusion  equations;  and  the  two-dimensional  steady- 
state  equation. 


AIR  POLLUTION  TRANSPORT 
MODELING 


I.  Introduction 

1.1  Background 

During  the  last  decade  there  has  been  an  incretising  public  interest  and  concern 
about  environmental  issues,  in  particular,  air  pollution.  Zannetti  (29)  stated  some  of 
the  issues  of  interest;  1)  the  “greenhouse”  effect,  which  could  cause  an  increase  in  the 
earth’s  average  temperature  due  to  increasing  concentrations  of  carbon  dioxide  in  the 
atmosphere;  2)  the  possible  depletion  of  the  ozone  layer,  a  natural  screen  of  harmful 
solar  radiation,  by  certain  pollutant  species  emitted  by  anthropogenic  activities;  3) 
indoor  air  pollution,  such  as  asbestos  and  radon  gas;  4)  nuclear  disasters,  especially 
after  the  Chernobyl  disaster;  5)  atmospheric  visibility  and  its  impairment  by  air 
pollutants;  and  6)  risk  assessment,  especially  in  the  prevention  of  accidental  releases 
of  toxic  pollutants.  One  of  the  simplest  air  pollution  transport  models  is  a  plume 
model  (more  specifically,  a  Gaussian  plume  model).  Many  air  pollution  plume  models 
have  been  developed  to  determine  the  concentrations  of  pollutants  in  the  air  and  the 
overall  quality  of  the  atmosphere.  These  plume  models  are  mostly  “specific”  in 
nature;  that  is  they  were  developed  with  a  specific  pollutant,  source  of  pollutant,  or 
change  of  concentration  in  mind. 

l.S  Problem 

Plume  models  have  been  used  extensively  in  determining  and  predicting  changes 
in  concentrations  of  air  pollutants.  Most  models  use  analytical  solutions  to  the  par¬ 
tial  differential  equations  used  in  pollution  modeling  in  their  calculations.  This  study 
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will  analyze  the  following  Gaussian  plume  models;  SCREEN  (10),  AFTOX  ( 19),  and 
GAUSPLUM  (the  progra-.  developed  in  this  research),  all  of  which  use  analytical 
solutions  in  the'r  cal'"  lations.  This  study  will  also  calculate  numerical  solutions  of 
the  advection,  steady  state,  and  advection-diffusion  equations  using  finite  difference 
methods  in  an  attempt  to  get  solutions  in  cases  of  low  wind  speed  which  the  Gaussian 
pi  11. lie  model  does  not  handle. 

1.3  Summary  of  Current  Knowledge 

Since  the  middle  1980s,  researchers  have  addressed  the  cause  and  effect  re¬ 
lationships  derived  from  air  pollution  models  in  general  (20).  The  results  of  this 
research  has  sparked  interest  in  studying  cause  and  effect  relationships  in  plume 
modeling  and  specifically  in  urban  areas.  Much  work  has  been  done  on  plume  rise 
and  diffusion  models,  (14,  15,  21)  and  on  modeling  pollution  in  an  urban  environ¬ 
ment  (9,  8).  There  have  also  been  studies  conducted  which  evaluated  and  verified 
mathematical  models,  and  these  studies  will  be  further  discussed  in  Chapter  11. 

1.4  Scope 

It  is  not  the  intent  of  this  study  to  completely  explain  all  plume  models  or  all 
air  pollution  transport  models.  Such  an  explanation  would  be  beyond  the  scope  of 
this  thesis.  The  thesis  is  limited  to  the  study  of  some  specific  plume  models.  The 
scope  of  this  research  is  limited  to  the  following: 

•  The  development  of  a  basic  pollutant  transport  model  using  the  Ada  program¬ 
ming  language. 

•  A  description  of  the  mathematical  a.spects  of  air  pollution  transport  models. 

•  An  analysis  and  comparison  of  the  numerical  results  of  the  models. 

•  Conclusion  reached  from  the  results,  and  any  recommendations  made  for  future 
research. 
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1.5  Approach 

The  basic  approach  taken  in  this  research  follows: 

1.  This  study  describes  the  general  mathematical  formula,  the  Gaussian  plume 
equation,  used  in  plume  models. 

2.  The  research  examines  the  specific  differences  among  the  models  and  how  the 
general  model  was  modified  a.s  a  result  of  such  differences.  Most  of  the  differ¬ 
ences  are  due  to  where  the  model  is  intended  to  be  used,  and  how  conservative 
the  user  wants  to  be.  Location,  terrain,  climate,  and  type  of  source  are  also 
criteria  specific  to  each  model. 

3.  Example  problems  will  then  be  applied  to  the  models  and  a  comparison  of  the 
results  will  be  conducted. 

4.  The  research  will  evaluate  several  finite  difference  schemes  and  apply  them  to 
solving  the  advection  equation,  the  one-,  two-,  and  three-dimensional  advec- 
tion-diffusion  equations,  and  the  two-dimensional  steady  state  equation. 

1.6  Materials  and  Equipment 

The  SCREEN  model  used  in  this  research  wzis  obtained  from  the  United  States 
Department  of  Commerce,  National  Technical  Information  Siervice,  Computer  Prod¬ 
ucts  Division  in  Springfield,  Virginia  via  the  United  States  Environmental  Protec¬ 
tion  Agency’s  Support  Center  for  Regulatory  Air  Models  Bulletin  Board  System 
(SCRAM-BBS)  with  the  Bulletin  Board  number  (919)  541-5742.  AFTOX  was  ob¬ 
tained  from  faculty  in  the  physics  department  at  the  Air  Force  Institute  of  Technology 
(AFIT). 

The  computer  facilities  at  AFIT  are  sufficient  to  run  the  software  models  used 
in  this  research.  No  other  computer  equipment  or  support  is  needed. 


1.1  Summary  of  Thesis 


The  organization  of  this  thesis  follows: 

Chapter  II  gives  an  overview  of  current  literature  concerning  plume  models, 
and  numerical  solutions  to  advection  and  diffusion  equations. 

Chapter  III  discusses  the  overall  methodology  of  this  ’•esearch  in  plume  mod¬ 
eling  including  the  numerical  methods  used  to  find  solutions  to  the  advection  and 
diffusion  equations. 

Chapter  IV  shows  the  results  of  the  research  and  discusses  the  comparisons 
made. 

Chapter  V  discusses  the  conclusions  derived  from  this  research  and  recommen¬ 
dations  for  mtuie  research. 

Appendix  A  contains  the  Ada  source  code  for  the  basic  plume  model,  GAUS- 
PLUM,  developed  in  this  research. 

Appendix  B  contains  the  Ada  source  code  for  the  advection  and  advection- 
diffusion  equations. 

Appendix  C  describes  the  one-dimensional  diffusion  equation,  and  the  three- 
dimensional  steady  state  equation. 
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II.  Literature  Review 


2. 1  Introduction 

This  literature  review  examines  current  research  in  the  areas  of  pollution  mod¬ 
eling  and  methods  of  finding  numerical  solutions  of  the  advection  and  diffusion  equa¬ 
tions.  Equation  (2.1)  shows  the  three  dimensional  advect ion-diffusion  equation  on 
which  most  models  are  based. 


dc  _dc  d‘^c  d'^c  d'^c 


(2.1) 


Various  approaches  to  air  pollution  plume  modeling  are  described  in  the  lit¬ 
erature  (14,  7).  This  review  addresses  the  Gaussian  and  Lagrangian  plume  models. 
The  Gaussian  plume  model  is  the  most  common  air  pollution  model  used  because 
of  its  relative  ease  of  use  with  easily  measurable  meteorological  parameters  (29). 
One  use  of  the  Gaussian  plume  equation  is  found  in  the  trajectory  plume  model. 
Okamoto’s  research  of  trajectory  plume  modeling  (22),  and  Schohl’s  research  of  La¬ 
grangian  plume  modeling  (24)  is  addressed  in  this  literature  review.  A  variety  of 
numerical  and  approximate  solutions  to  plume  modeling  are  also  described  in  the 
literature  (13,  17,  16,  23).  This  literature  review  will  also  look  at  some  numerical 
methods  and  some  of  the  limits  of  pollution  modeling. 

Two  terms  that  need  further  definition  are  stability  condition  and  stability 

class: 

A  stability  condition  is  a  constraint  that  is  put  on  variables  in  a  solution 
method.  When  the  value  of  the  variables  are  in  the  range  of  the  constraint  then 
small  changes  in  the  variables  produces  small  changes  in  the  results  and  the  method 
is  said  to  be  stable  (3). 

A  stability  class  is  a  parameter  used  to  represent  the  atmospheric  conditions 
(ranging  from  extremely  unstable  to  extremely  stable,  A  to  G).  It  is  determined 
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using  the  standard  deviation  of  the  horizontal  wind  direction  fluctuation  and  the 
wind  speed  at  10  meters  in  meters  per  second.  See  Turner’s  workbook  (27)  or  Table 
3  in  the  AFTOX  User’s  Guide  (19)  for  classifications. 

Stability  condition  and  stability  class  are  common  terms  used  in  the  analysis 
of  air  pollution  transport  modeling. 

2.2  Trajectory  plume  model 

Among  the  different  models  used  for  calculating  pollutant  concentrations,  a 
Gaussian  plume  model  is  the  most  common.  The  Gaussian  plume  model  is  based 
on  a  formula  that  describes  the  three-dimensional  pollutant  concentration  generated 
by  a  point  source  under  stationary  meteorological  and  emission  conditions  (22,  27). 
Its  calculation  error  becomes  a  serious  problem  under  weak,  variable-direction  wind 
conditions.  It  is  for  this  reason  that  Okamoto’s  (22)  research  of  a  trajectory  plume 
model  or  a  plume  segment  model  was  conducted.  This  model  treats  time  varying 
transport  conditions  and  changes  in  wind  direction  and  speeds.  In  the  segmented 
plume  approach,  the  plume  is  broken  up  into  independent  elements  (plume  segments 
or  sections)  whose  initial  features  and  time  dynamics  axe  a  function  of  emission 
conditions  and  local  meteorological  conditions  encountered  by  the  plume  segment 
(29:165-7).  Segments  are  sections  of  a  Gaussian  plume.  Therefore,  the  concentra¬ 
tion  of  air  pollution  for  each  segment  can  be  calculated  using  the  Gaussian  plume 
model.  The  trajectory  plume  model  then  adds  segments  together  to  get  concentra¬ 
tion  distributions  over  the  entire  plume. 

2.3  Lagrangian  plume  model 

Another  model  which  uses  segmented  plumes  to  calculate  concentrations  of  air 
pollutants  is  the  Lagrangian  box  model.  The  Lagrangian  box  model  breaks  the  plume 
into  boxes  (segments),  each  of  which  follows  the  average  wind  flow.  Because  the  boxes 
move  with  the  average  wind  flow,  the  Lagrangian  box  model  provides  concentration 
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outputs  along  trajectories  and  consequently  comparisons  with  concentrations  at  fixed 
grid  points  of  traditional  finite  difference  or  finite  element  solutions  are  difficult  (29). 
Schohl  (24)  discusses  an  interpolation  method  to  use  in  a  Lagrangian  scheme: 

“Lagrangian  schemes  for  numerical  approximation  of  concentration  a- 
mounts  are  capable  of  providing  satisfactory  accuracy  with  minimal  com¬ 
putational  effort.  The  overall  accuracy  of  these  schemes  is  dependent  on 
the  interpolation  method  used  to  obtain  the  approximation.”  (24) 

Schohl’s  (24)  research  shows  that  the  cubic-spline  interpolation  method  is  one  which 
provides  better  than  average  accuracy  to  the  Lagrangian  model.  For  more  informa¬ 
tion  on  cubic  spline  interpolation  see  references  (3:126-7)  and  (28).  This  interpolation 
provides  a  fourth  order  approximation  with  minimal  computational  effort. 


2-4  Numerical  and  approximate  solutions 

2.4-1  Plume  rise  model  solution  method.  Krishnamurthy  (18)  states  that 
little  seems  to  have  been  published  which  gives  useful  guidelines  concerning  the 
accuracy  of  approximations  used  in  air  pollution  transport  models.  In  fact,  he  says 

“In  general,  the  application  of  plume  models  requires  numerical  solu¬ 
tions  of  the  governing  conservation  equations  which  are  highly  non-linear. 
Rather  than  obtain  such  ’’exact”  solutions  it  is  easier  to  use  various  ap¬ 
proximations  to  them.  Sometimes,  however,  the  accuracy  of  such  ap¬ 
proximations  relative  to  the  exact  solution  is  quite  uncertain.”  (18:2083) 

Krishnamurthy’s  research  calculates  numerical  solutions  using  the  plume  rise  model 
of  Hoult,  Fay,  and  Forney  (15).  Various  eisymptotic  approximations  are  assessed  over 
a  wide  range  of  parameters.  Mass,  momentum,  energy,  turbulence,  temperature, 
and  wind  speed  and  direction  are  the  parameters  used  in  these  approximations. 
Comparisons  of  the  numerical  and  approximate  solutions  are  made  and  shown  to 
agree  over  a  fairly  wide  range  of  the  parameters. 

2.4.2  Advection  equation  solution  methods.  Many  methods  are  described 
in  the  literature  which  numerically  solve  partial  differential  equations.  In  partic- 
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ular,  Chock,  in  “A  Comparison  of  Numerical  Methods  for  Solving  the  Advection 
Equation-I,  II,  and  III”  (4,  5,  6)  compares  several  algorithms  and  variations  of  each 
which  are  used  for  solving  the  two-dimensional  advection  equation  which  is 


dc  dc  dc 


(2.2) 


where  c  is  the  concentration,  and  u  and  v  are  the  x  and  y  components  of  the  wind 
speed. 


The  methods  are  compared  in  terms  of  accuracy,  speed,  and  storage  require¬ 
ments.  Recommendations  are  made  in  the  articles  as  to  how  useful  the  methods 
would  be  for  air  quality  modeling. 


2.4-3  Diffusion  equation  solution  methods.  There  are  also  many  techniques 
in  the  literature  for  solving  th  itmospheric  diffusion  equation 


dc  d'^c  .  d^c 

where  in  equation  (  2.3)  c  is  the  concentration  of  the  pollutant  in  question  and 
kx,  ky,  and  k^  are  the  x,y,  and  z  diffusivity  terms  respectively.  McRae,  Goodin,  and 
Seinfeld  look  at  several  techniques  for  solving  both  the  advection  equation  (2.2)  and 
i  ' :  diffusion  equation  (2.3)  in  their  one-dimensional  form  in  “Numerical  Solution  of 
the  Atmospheric  Diffusion  Equation  for  Chemically  Reacting  Flows”  (21). 

2.4.4  Advection-Diffusion  solution  methods.  The  one-dimensional  advec- 
tion-diffusion  equation,  or  as  Strikwerda  (26:129-31)  calls  it,  the  convection-  diffusion 
equation,  is 


dc  _dc  d^c 


(2.4) 
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where  c  is  the  concentration,  u  is  the  advection  coefficient  and  kx  is  the  difFusivity 
coefficient.  Strikwerda  looks  at  this  equation  two  ways;  one  way  using  a  substitution 
oiy  =  X  —  ut  with  w{t,y)  =  c{t,y  +  ut),  and  the  other  way  using  the  forward-time 
central  space  finite  difference  method  (see  Section  3.4.2  for  this  method). 

2.5  Limits  of  air  pollution  modeling 

Numerical  and  approximate  solutions  may  agree  over  a  wide  range  of  parame¬ 
ters,  but  the  limits  of  meteorological  modeling  can  lead  to  uncertainty  in  atmospheric 
parameters  such  as  turbulence,  wind  and  temperature.  This  can  lead  to  uncertainty 
in  the  models  that  use  these  parameters.  The  Gaussian  plume  model  doesn’t  per¬ 
form  well  when  the  wind  is  weak,  and  an  uncertainty  in  the  wind  parameter  can 
add  to  this  poor  performance.  Benarie  (1)  states  that  “turbulence  is  a  meteorolog¬ 
ical  quantity  that  can  only  be  approximated  and  then  only  in  the  most  ideal  cir¬ 
cumstances”.  Meteorological  conditions  are  inherently  variable  and  this  variability 
causes  an  uncertainty  in  the  accuracy  of  meteorological  parameters.  Plume  models 
which  incorporate  these  parameters  inherit  this  uncertainty.  Therefore,  numerical 
and  approximate  solutions  of  air  pollution  models  are  dependent  on  the  limits  of 
meteorological  modeling  used  to  calculate  the  parameters  for  pollution  models. 

2.6  Conclusion 

Trajectory  plume  i  o  leling  and  Lagrangian  plume  modeling  both  use  segmen¬ 
tation  of  pollution  plumes  in  their  calculation  of  pollutant  concentrations.  The 
trajectory  model  uses  the  Gaussian  plume  equation  to  calculate  the  concentra¬ 
tion  in  each  segment.  The  Lagrangian  model  uses  interpolation  methods  for  cal¬ 
culations  between  segments.  Many  methods  for  solving  advection,  diffusion,  and 
advection-diffusion  equations  used  in  these  models  are  found  in  the  literature  includ¬ 
ing,  the  plume  rise  model,  several  advection  equation  solution  methods  in  Chock’s 
research,  techniques  to  solve  both  oiK-  .L.nensional  advection  and  diffusion  equations 
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in  McRzie’s  research,  and  methods  used  by  Strikwerda  to  solve  the  one- dimensional 
advection-difFusion  equation.  However,  the  accuracy  of  these  numerical  and  approx¬ 
imate  solutions  are  dependent  on  the  parameters  used  in  the  models.  Much  more 
research  is  needed  in  the  areas  of  numerical  and  approximate  solutions  and  the  de¬ 
termination  of  the  parameters  used  in  pollution  models.  This  thesis  will  use  finite 
difference  methods  to  solve  these  equations. 
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III.  Methodology 


3. 1  Introduction 

The  purpose  of  this  study  is  to  examine  air  pollution  transport  models  and  the 
partial  differential  equations  used  in  air  pollution  modeling.  This  research  examines 
the  following  models:  SCREEN  (10),  AFTOX  (19),  the  program  GAUSPLUM  de¬ 
veloped  in  this  study,  and  numerical  schemes  used  to  solve  the  advection,  diffusion, 
two-dimensional  steady  state,  and  one,  two  and  three  dimensional  advection-diffusion 
equations. 

The  programming  language,  a  brief  description  of  the  capabilities,  and  the 
mathematical  equations  used  in  SCREEN,  AFTOX,  and  GAUSPLUM  tire  discussed. 
The  two  primary  capabilities  used  in  this  research  include  finding  the  concentration 
of  a  pollutant  at  a  given  location,  and  finding  the  location  and  value  of  the  maximum 
concentration.  A  comparison  of  these  capabilities  in  each  model  is  done  using  six 
example  problems  taken  from  the  Workbook  of  Atmospheric  Dispersion  Estimates 
(27).  These  problems  are  described  in  detail  below. 

The  Lax-Friedrichs  and  leapfrog  finite  difference  methods  will  be  used  on  the 
one-dimensional  advection  equation.  The  forward-time  central-space  finite  difference 
method  will  be  used  on  the  advection-diffusion  equations. 

3.3  Description  of  Gaussian  plume  models 

This  section  describes  the  SCREEN  and  AFTOX  models,  and  the  GAUSPLUM 
model  developed  in  this  research. 

3.2.1  The  SCREEN  model.  SCREEN  is  written  in  FORTRAN  program¬ 
ming  language.  SCREEN  estimates  pollutant  concentration  from  continuous  sources 
using  a  Gaussian  plume  model  that  incorporates  source-related  factors,  such  as  emis¬ 
sion  rate,  stack  gas  temperature,  stack  height,  stack  inside  diameter,  and  stack  gas 
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exit  velocity,  and  meteorological  factors,  such  as  ambient  temperature,  wind  speed 
and  direction.  The  Gaussian  model  equations  are  described  in  Turner’s  workbook 
(27). 


The  basic  equation,  for  determining  ground-level  concentrations  under  the  cen¬ 
terline  of  the  plume,  used  in  the  SCREEN  model  (10)  is; 


c  =  [9/(27rUj<Ty<T^)]  •(  exp[-i((0r  -  K)la^f] 

-I-  exp[-|((2:,.  -1-  he)la^f] 

+  [  exp[-i((2r  -h,-2N z^)/(r^)^] 

+  exp[-i((2r  +  he-  2Nzi)/a^y] 
-|-exp[-i{(2,.  -  he  +  2Nzi)/a^y] 

-H  exp[-|((2,.  +  he  +  2Nzi)/a^y]]) 


where 

c  =  concentration  {gfm^) 
q  =  emission  rate  (g/s) 

TT  =  3.14159 

Uj  =  stack  height  wind  speed  (m/s) 

<7y  =  lateral  dispersion  parameter  (m) 

<72  =  vertical  dispersion  parameter  (m) 

Zr  =  receptor  height  above  ground  (m) 
he  =  plume  centerline  height  (m) 

Zi  =  mixing  height  (m) 

k  =  summation  limit  for  multiple  reflections  of  plume 
off  of  the  ground  and  elevated  inversion,  usually  <  4 

This  equation  accounts  for  the  multiple  eddy  reflections  from  both  the  ground 
and  the  stable  layer  (2^)  and  was  suggested  by  Bierly  and  Hewson  (2).  The  derivation 
of  equation  (3.1)  is  discussed  in  Section  3.2.3. 

3.2.2  The  AFTOX  model.  AFTOX  is  a  program  written  in  the  Basic  pro¬ 
gramming  language.  The  two  parts  of  the  USAF  Toxic  Chemical  Dispersion  Model 
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(AFTOX)  used  in  this  research  are  the  calculation  of  the  toxic  chemical  concentration 
at  a  specific  location  aria  the  strength  and  location  of  the  maximum  concentration. 
Refer  to  AFGL-TR-88-0009  (19)  for  further  uses  of  the  model.  AFTOX  uses  the 
Gaussian  puff  equation,  and  the  Gaussian  plume  equation  in  its  calculations.  These 
equations  2issume  a  Gaussian  distribution  of  concentration  and  conservation  of  the 
pollutant  during  transport  and  diffusion. 

The  AFTOX  model  uses  three  basic  models,  the  Gaussian  puff  equation  (see 
eq.3.2),  the  Gaussian  puff  equation  when  an  inversion  is  present  (see  eq.3.3),  and  the 
Gaussian  plume  model  (see  eq.3.4). 

The  Gaussian  puff  model  is 

c{x,  y,z,t-i)=  \q{i)j{{‘lTrfna^aya^)] 

•exp[-i((ar-u(<-f))/<y^)2] 

•exp[-(j^/<7j,)V2] 

•(exp[-|((2  -  H)la^f]  +  exp[-i((2  +  H)la^f]) 


where 

c  is  concentration  in  the  puff  at  (x,y,z)  at  {t  —  i) 

q  is  total  mass  of  the  puff 

u  is  wind  speed  at  10m 

(Tx  is  downwind  dispersion  parameter 

cTy  is  lateral  dispersion  parameter 

(Tz  is  vertical  dispersion  parameter 

t  is  total  elapsed  time  of  pollution  emission 

t  is  time  of  puff  emission 

{t  —  t)  is  elapsed  time  since  puff  emission 

H  is  height  of  the  source. 
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The  Caussian  puff  model  with  an  inversion  has  the  addition  of  the  following 
expression  to  the  last  two  terms  of  equation  (3.2) 

E/v(  expl-l((2-//-2Wi)/<T.)^|+expl-l((2  +  //-2Ari)/,T,)"l+ 

+  2NL)la^f]  +  exp[-i((x  +  //  +  2NL)ja^f\). 

where  L  is  the  mixing  layer  height  and  N  is  the  number  of  reflections  caused  by  the 
inversion.  This  equation  is  similar  to  equation  (3.1)  used  in  the  SCREEN  model. 

The  Gaussian  plume  model  used  in  AFTOX  is 


c  = 


1.  y 


\,z-H 


2ir<Tx(TyU 


1  2  + 


exp[-5(^)^]  -  (exp[-l(:i^)^]  -f  exp[-^(:^)-^]).  (3.4) 


2'  a 


2'  a. 


The  derivation  of  equations  (3.2),  (3.3),  and  (3.4)  is  discussed  in  Section  3.2.3. 


3.2.3  The  GAUSPLUM  model  from  this  research.  GAUSPLUM  is  written 
in  the  Ada  programming  language.  It  uses  the  Gaussian  plume  equation  with  the 
Pasquill-GifiFord  (Ty  and  cr^  described  by  Zannetti  (29:149-50). 

The  basic  Gaussian  plume  equation  used  in  GAUSPLUM  is: 


c  = 


2TC<Ty<7gU 


^  \2l  „„„r  ^  \2 


exp[-x(— )  ]exp[--(- 

l  <7, 


(3.5) 


where  c  is  the  concentration  at  r  =  {xr^yri^r)  due  to  emissions  from  the  source 
at  (x,,  j/,,  2,);  q  is  the  emission  rate;  ay  and  Oz  are  the  horizontal  and  vertical, 
respectively,  dispersion  parameters;  u  is  the  horizontal  wind  speed;  and  hg  is  the 
effective  emission  height.  The  coordinate  system  in  this  model  has  the  x-axis  in  the 
u  direction. 


Equations  (3.1),  (3.2),  (3.3),  (3.4),  and  (3.5)  can  be  derived  from  the  three- 
dimensional  advection-diffusion  equation  which  is 
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(3.6) 


“1“  UCx  —  ICx^xx  ^y^yy 

In  equation  (3.6),  at  steady  state,  cj  =  0  so  the  equation  reduces  to 


tlCx  —  ^x^xx  d"  ^y^yy  d“  ^z^zz*  (3**^) 

For  many  air  pollution  transport  problems,  the  kxCxx  term  is  negligable  compared 
to  the  UCx  (25:542-543).  Then  equation  (3.7)  reduces  to 


UCx  -  kyCyy  “I”  Jc^C^x*  (3*3) 

If  u,  ky,  and  k^  are  constant  and  the  source  is  a  point  source,  an  exact  closed 
form  solution  of  equation  (3.8)  can  be  obtained  using  Fourier  transform  techniques 
(25:556).  This  solution  is 


S'- fc*'- 


(3.9) 


See  Seinfeld  page  543  and  556  for  details.  If  (<7^u)/(2i)  and  (<7^u)/(2x)  are  both  con¬ 
stant,  then  letting  ky  =  (a-y«)/(2x)  k^  =  (<7^u)/(2x)  and  substituting  into  equation 
(3.9)  gives 


c  = 


2TrCTy(TxU 


2V, 


(3.10) 


Equation  (3.10)  is  the  same  as  equations  (3.1),  (3.2),  (3.3),  (3.4),  and  (3.5)  except 
they  include  a  term,  H  or  /le,  in  the  z  term  to  account  for  the  height  of  the  source. 
Equation  (3.1)  also  includes  a  term  for  when  an  inversion  is  present  in  the  atmosphere 
and  cissumes  y  =  0  so  there  is  no  y  term  in  the  equation.  It  is  these  equations  which 
can  be  derived  from  the  three-dimensional  advection-dilfusion  equation  that  are  used 
in  the  SCREEN,  AFTOX,  and  GAUSPLUM  models  to  solve  the  following  problems. 
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3.S  Description  of  example  problems 

This  section  describes  the  problems  taken  from  Turners  workbook  used  to 
compare  the  three  models. 

•  Problem  1  from  the  workbook  shows  a  ground-level  calculation  directly  down¬ 
wind  at  a  distance  of  3000  meters.  The  ground-level  source  emits  3  gjsec  of 
oxides  of  nitrogen  with  no  effective  rise.  It  is  an  overcast  night  with  a  7  mfsec 
wind  speed.  This  indicates  a  stability  class  D. 

•  Problem  2  from  the  workbook  calculates  a  concentration  500  meters  directly 
downwind  from  a  source  with  an  effective  height  of  60  meters.  It  is  an  overcast 
winter  morning  at  0800.  The  source  emits  an  estimated  80  g/sec  of  sulfur 
dioxide  into  a  wind  of  6  mfsec.  The  stability  class  is  still  D. 

•  Problem  3  has  the  same  conditions  as  problem  2  at  the  same  distance  downwind 
but  at  a  distance  of  50  meters  from  the  x-axis.  The  SCREEN  program  only 
calculates  concentrations  directly  downwind  from  the  source  so  it  does  not 
apply  to  this  problem. 

•  Problem  4  has  a  source  emitting  151  y/sec  at  an  effective  height  of  150  meters. 
This  problem  asks  for  the  distance  and  the  value  of  the  maximum  ground-level 
concentration  on  a  sunny  summer  afternoon  with  a  10  meter  wind  speed  of  4 
mfsec  from  the  northeast.  This  is  a  class-B  stability. 

•  Problem  5  has  the  same  conditions  as  problem  4  except  it  is  on  an  overc^lst 
day.  The  stability  class  becomes  D  for  this  problem. 

•  Problem  11  in  the  workbook  has  the  same  conditions  as  problem  4  except  it 
asks  for  the  distance  and  value  of  the  maximum  concentration  on  the  plume 
centerline  on  a  clear  night  with  a  wind  speed  of  4  mfsec.  This  give  stability 
class  E. 
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These  problems  will  be  used  to  compare  SCREEN,  AFTOX,  and  GAUSPLUM 
while  finite  difference  schemes  will  be  used  to  solve  the  following  mathematical  mod¬ 
els  and  compare  the  solutions  to  exact  analytical  solutions. 

3.4  Mathematical  models 

This  section  describes  the  mathematical  finite  difference  schemes  examined  in 
this  research.  The  advection  equation,  the  one,  two  and  three  dimensional  advection- 
diffusion  equations,  and  the  two  and  three  dimensional  steady-state  equations  are 
described  here. 

3.4.1  The  advection  equation:  The  advection  equation  (see  eq.  3.11),  also 
known  as  the  one-way  wave  equation,  is  a  hyperbolic  partial  differential  equation. 
This  subsection  describes  the  simple  advection  equation  and  the  numerical  methods 
used  to  solve  it  in  this  research.  The  advection  equation  is 


dc  _dc 


(3.11) 


also  written  as  q  -f  ucx  =  0,  where  the  subscript  denotes  differentiation,  i.e.,  C(  = 
dcjdt,  and  u  is  the  average  wind  speed. 

Two  finite  difference  schemes  are  used  on  the  following  initial  boundary  value 
problem: 


Ct  +  Cx  =  0  on  —  2  <  X  <  3, 0  <  < 


with  initial  data 


Co(a:) 


1  —  |x|  if  )x|  <  1 
0  if  |x]  ^  1 


The  boundary  condition  at  time  t  is 


(3.12) 


(3.13) 
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c{x,t)  =  0 


(3.14) 


when  X  =  —2,  or  x  =  3. 


The  finite  difference  schemes  to  be  used  on  the  above  initial-boundary  value 
problem  are  the  Lax-Friedrichs  scheme  (3.15)  and  the  leapfrog  scheme  (3.16)  as  found 
in  Strikwerda  (26:13-4).  The  notation  is  the  same  as  c(<„,x,n)  or  in  other  words 
it  is  the  value  of  c  at  the  grid  point  (<„,Xm). 


+  C-i)  , 


At 


+  w- 


m+l 


—  c: 


m— 1 


2Ax 


=  0 


(3.15) 


-n+1 


—  c: 


,n—l 


2At 


^Tl 

+  u^ 


—  c: 


m— 1 


2Ax 


=  0 


(3.16) 


where  u  =  1  for  this  problem,  At  is  the  time  incrementer,  Ax  is  the  space  incre- 
menter,  m  and  n  are  integer  grid  counters  for  x  and  t  respectively.  Solving  the 
schemes  for  gives  a  linear  combination  of  c  at  levels  n  and  n  —  1.  The  Lax- 
Friedrichs  scheme  (3.15)  can  be  written  eis 


(3.17) 


where  A  =  At/ Ax.  Likewise,  the  leapfrog  scheme  (3.16)  can  be  written  as 


C'  =  C  -  «A(c” +1  -  c”  .i)  (3.18) 

again  with  A  =  At/Ax,  and  u  =  1. 

The  stability  condition  for  the  Lax-Friedrichs  method  is  uA  <  1,  or  since  u  =  1 
in  this  case,  A  <  1.  It  is  found  by  replacing  c”  with  in  equation(3.17)  and 

solving  for  g  which  yields 
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g  =  (e-®  +  e-‘‘')/2  -  «A(e’®  -  e-‘®)/2  (3.19) 

where  i  =  (—1)'^^.  The  quantity  g  is  called  the  amplification  factor.  This  is  equiva¬ 
lent  to 

g  =  cosO  —  iuXsmO.  (3.20) 

The  stability  condition  |gi|  <  1  comes  from  Theorem  2.2.1  in  Strikwerda  (26:42). 
Note  that  Ipp  is 


kP  =  cos^  ^ -f  sin^  (3.21) 

Therefore,  |5f|  <  1  if  |uA|  <  1.  Thus  the  stability  condition  uA  <  1  since  both  it  and 
A  are  positive  in  this  case. 

If  u  =  1,  and  A  =  0.8  the  stability  condition  is  satisfied  for  both  schemes. 
Consistency  occurs  when  the  local  truncation  error  vanishes  as  Ax  0  and  At  -+  0 
(28:606).  The  local  truncation  error  is  the  difference  between  the  solution  of  the 
difference  equation  at  a  point,  and  the  solution  of  the  differential  equation  at  the 
same  point.  By  substituting  Taylor  series  expansions  into  the  Lax- Friedrichs  and 
leapfrog  schemes  it  can  be  shown  that  both  schemes  are  consistent  (26:21-23).  Both 
conditions,  stability  and  consistency,  must  be  met  for  the  scheme  to  be  convergent. 
Both  the  Lax-Friedrichs  scheme  and  the  leapfrog  scheme  are,  therefore,  convergent 
for  A  =  0.8  since  both  stability  and  consistency  conditions  are  met.  Figures  (4.3) 
and  (4.4)  in  Section  4.3.1  show  that  both  schemes  are  convergent  for  A  =  0.8. 

When  A  =  1.6  is  used  for  the  Lax-Friedrichs  method  it  is  not  convergent.  This 
is  because  the  stability  condition  is  not  met.  However,  as  Ax  — v  0  and  A<  — >  0 
the  solution  of  the  scheme  does  approach  the  solution  of  the  advection  equation,  as 
long  as  >  0,  so  the  scheme  is  consistent  (see  Strikwerda,  example  1.4.2,  page 
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21.)(26:21).  The  results  of  using  these  finite  difference  schemes  will  be  discussed  in 
Section  4.3.1. 

3.4.2  The  l-D  advection- diffusion  equation:  The  one- dimensional  advec- 
tion-diffusion  equation  (3.22)  is  also  known  as  the  one-dimensional  convection-diffu¬ 
sion  equation, 


C(  -}-  UCx  —  ^x^xx  (3.22) 

where  kx  is  a  positive  number  and  we  assume  that  u  is  also  positive.  The  forward¬ 
time  central-space  scheme  used  to  solve  equation  (3.22)  is 


_ 

At 


+  u 


"m+l 


—  C 


m— 1  _ 


2Ax 


=  kx 


jn+i 


-  2cZ  +  c] 


m— 1 


(3.23) 


which  is  equivalent  to 


=  (1  -  2kxn)cl^  +  k^ti{l  -  a)C+i  -|-  kxfi{l  +  a)c”  _i  (3.24) 

where  n  =  Atj  Ax^  and  a  —  uAz/2fcj,.  The  finite  difference  scheme  is  used  on  the 
following  initial-boundary  value  problem: 


Cf  -|-  UCx  —  kxCxx 


on  —  2<z<3,0<< 


with  initial  data 


Co(a:) 


1  —  |x|  if  |z|  <  1 
0  if  |z|  >  1. 


The  boundary  condition  at  time  t  is 


(3.25) 


(3.26) 


c(z,  t)  =  0 


(3.27) 
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when  X  =  —2,  or  a;  =  3. 

The  stability  condition  for  this  method  is  kxH  <  1/2.  The  stability  analysis 
is  similar  to  that  done  on  the  advection  equation  above,  i.e.  by  replacing  c”  with 
gTi^xme  jjj  equation  (3.23)  and  solving  for  g  gives 

g  =  I  —  4kxfi  sin^  —  iuX  sin  6.  (3.28) 

Using  the  condition  <  1  from  Theorem  2.2.1  and  Theorem  2.2.3  and  Corollary 
2.2.2  in  Strikwerda  (26),  which  shows  that  the  first  derivative  term  can  be  ignored, 
gives  ikxfism^  <  2  and  so  kxfi  <  1/2.  Thus  this  method  is  conditionally  stable. 
The  results  of  using  this  finite  difference  scheme  to  solve  equation  (3.22)  will  be 
presented  in  Section  4.3.2. 

3.4. S  The  2-D  advection-diffusion  equation:  The  two-dimensional  advec- 

tion-diffusion  equation  (see  eq.  3.29)  includes  diffusion  in  both  the  downwind  (x- 
axis)  ajid  cross  wind  (y-axis)  directions.  Advection  is  assumed  to  be  negligable  in  the 
crosswind  direction  since  the  wind  direction  is  in  the  x-axis  direction.  Therefore,  the 
only  advection  term  is  ucx- 


Cf  “1“  U/Cx  —  kxCxx  d”  (3.29) 

where  u  is  the  wind  speed  (positive  constant),  and  kx  =  ky  are  constant  in  this 
research  and  are  the  x-axis  and  y-axis  diffusivity  terms  respectively.  The  forward¬ 
time  central-space  finite  difference  scheme  used  to  solve  equation  (3.29)  is 


^n+l  _  „n 

At 


2Ax 


{kx/Ax  )(c^4i_p  p  -|-  cJ)j_j  p) 

{ky/Ay^)(cl^y^,  -  2c”  ,p  +  c;j,,p_,)  (3.30) 
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which  is  equivalent  to 


Cmtp  =  C,p  -  ^“'^(C+I.P  -  C-l,p) 

+  *xa(C+i,p  -  2C,p  +  C"  _i_p) 

+  M(C.p+i-2C,p  +  c;;.,p_,)  (3.31) 

where  cJJ^p  =  c(xm,yp,i„)  =  c{TnAx,pAy,nAt),  X  =  At /Ax,  a  =  At/Ax^  and 
^  =  A</Aj/^ 

The  finite  difference  scheme  is  used  on  the  following  initial-  boundary  value 
problem: 


ct  +  ucx  =  kxCxx  +  kyCyy  ou  — 2<x<3,  —2<y<2,  0  <t  (3.32) 


with  initial  data 


CQix,y) 


< 


1-R  if  0  <  X  <  1 
0  otherwise 


where  R  =  sqrt(x*  -t-  y^).  The  boundary  condition  at  time  t  is 


(3.33) 


c{x,y,t)  =  0  (3.34) 

when  X  =  —2,  or  x  =  3,  or  y  =  ±2. 

The  stability  condition  for  this  method  is  kxfi  <  1/4.  The  stability  analysis  is 
similar  to  that  done  on  the  one-  dimensional  equation  above,  i.e.  by  replacing  c” 
with  in  equation  (3.30)  and  solving  for  g  gives 


g  =  I  —  ik 


Ak, 


:yfis\n'^^0 


■  iii.X  sin  9 


(3.35) 
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or  since  kx  =  ky 


g  =  I  —  Skxfisin^ -0  —  iuXsinO.  (3.36) 

£t 

Using  the  condition  1^1  <  1  from  Theorem  2.2.1,  Theorem  2.2.3  and  Corollary  2.2.2 
in  Strikwerda  (26),  which  shows  that  the  first  derivative  term  can  be  ignored,  gives 
Skxfi  <  2  so  kxH  £  1/4.  Thus  this  method  is  conditionally  stable.  The  results 

of  using  this  finite  difference  scheme  to  solve  equation  (3.29)  will  be  presented  in 
Section  4.3.3. 

3.4-4  The  3-D  advection-diffusion  equation:  The  three-dimensional  advec- 
tion-diffusion  equation  (3.37)  includes  diffusion  in  the  downwind  (x-axis),  crosswind 
(y-axis)  and  vertical  (z-axis)  directions.  Advection  is  assumed  to  be  negligable  in  the 
crosswind  and  vertical  directions  since  the  wind  direction  is  in  the  x-axis  direction. 
So  again,  the  only  advection  term  is  tiCx, 


“1“  uCx  —  kxOxx  d”  kyCyy  -f-  k^c^z  (3.37) 

where  u  is  the  wind  speed  (positive  constant),  and  kx  =  ky  =  kz  are  constant  in  this 
research  and  are  the  x-axis,  y-axis  and  z-axis  diffusivity  terms  respectively.  In  this 
research  two  variations  of  three-dimensional  equation  are  examined.  One  variation 
uses  a  forward-time  central-space  finite  difference  scheme  to  solve  equation  (3.37) 
with  the  condition  at  t  =  0 

c(i,i/,2, 0)  =  80  for  x  =  j/  =  z  —  0  (3.38) 

where  80  was  chosen  from  example  problem  2  described  above.  The  boundary  con¬ 
dition  at  time  t  is 
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c{x,y,z,t)  =  0 


(3.39) 


when  X  =  ±10  or  y  =  ±10,  or  z  =  ±10.  The  forwaxd-time  central-space  finite 
difference  scheme  used  to  solve  equation  (3.37)  is 


.n+l 


—  C 


i,j,k 


At 


±  u 


2Ax 


+  (*./A2’)(c?^jj,.,  -  2c’j  i,  +  cf  j  i., 


) 

)(3.40) 


which  is  equivalent  to 


^?,j,k  l/2t*A(c"^ij,j|,  —  c"_j  j|  jt) 

+  f^M<+i,j,k  -  Kj,k  +  fc) 

±  ^v^(c"j^l  —  2c” J  ±  ^?,j-l,k) 

+  ^Mcli,k+i  -  '^<j,k  +  <^lj,k-i) 


(3.41) 


where  ^  =  c{xi,yj,  Zk,tn)  =  c{iAx,jAy,kAz,nAt),  X  =  At/ Ax,  a  =  At/Ax^, 
(3  =  At/Ay^,  and  7  =  At/Az^. 

The  solution  to  this  scheme  will  be  compared  to  the  exact  analytical  solution 
from  Seinfeld  (25:536)  which  is 


c{x,  y,  z,  t)  =  [s/{S{irtf^^{kxkyk^y^'^)] 

•exp[-(x  -  uty/{Akj)  -  y‘^/{4kyt)  -  zV(4^zt)] 


(3.42) 


where  s  is  the  source  term,  or  80  in  this  case. 
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The  other  variation  uses  a  forward-time  central-space  finite  difference  scheme 
to  solve  equation  (3.43)  which  adds  a  source  term,  or  forcing  function,  ^(f)  •  /(i,  y,  z) 
to  equation  (3.37)  instead  of  initial  conditions  which  is 


Ct  UCx  —  kxCxx  d"  ^y^yy  d"  ^z^zz  d”  ^(f)  *  J  J/j  ^)'  (3.43) 


In  equation  (3.43) 


f{x,y,z)  =  sin 


{x  d-  10)7r  .  {y  -|-  10)7r  .  {z  -f  10)7r 


20 


sin 


20 


sin 


20 


(3.44) 


when  X  =  y  =  0,  and  z  =  height,  where  height  is  the  height  of  the  source,  and  6{t) 
is  the  delta  function. 

The  forward-time  central-space  finite  difference  scheme  used  to  solve  equation 
(3.43)  is 


^n+l 


-  C 


At 


—  c 


2Ax 


{kx! Ax  )(c"^i,j,jt  -|-  c”_i  j  j.) 

d-  {kylAy'^){clj^^^k  -  2c”j,*  d- 

d-  (A:^/Az*)(c^j-fc+i  -  2cl^  k  +  <i,*-i) 

+  6{t)-  f{x,y,z)  (3.45) 


which  is  equivalent  to 


-  l/2wA(c,\i,^.fc 


—  c, 


i-l,j,k) 


+  kxa{c^+ij,k  -  d- 
d-  ^y^(c"j+l,fc  ~  2c"^;t  d- 
+  kzl{clj,k+i  -  2<j,fc  d- 


) 

) 

) 
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+  8{t)-  f{x,y,z)-  At 


(3.46) 


where  c”^  =  c{xi,yj,  Zk^tn)  =  c{iAx,jAy,kAz,nAt),  A  =  At/ Ax,  a  =  At/Ax^, 
/3  =  At/Ay^,  7  =  At/Az^,  and  f{x,y,z)  is  the  forcing  function  or  source  term. 

The  solution  to  this  scheme  will  be  compared  to  the  exact  analytical  solution 
which  is 


c{x,  y,  z,  t)  =  exp[cj<]  •  sin[7r(a:  +  +  10)/20] 

•  sin[7r(j/  +  l0)/20]  •  sin[7r(2  +  10)/20] 

where  u  =  —{ir/i00){kx  +  ky  +  k^). 


(3.47) 


The  stability  condition  for  these  methods  is  kxfi  <  1/8,  when  kx  =  ky  =  k^. 
The  stability  analysis  is  similar  to  that  done  on  the  two-dimensional  equation  above, 
with 


g  =  I  —  16kxfisin^  —  ia\sin9.  (3.48) 

Using  the  condition  |^|  <  1  from  Theorem  2.2.1,  Theorem  2.2.3  and  Corollary  2.2.2 
in  Strikwerda  (26),  which  shows  that  the  first  derivative  term  can  be  ignored,  gives 
16^1^  sin^  and  kxH  <1/8.  Thus  this  method  is  also  conditionally  stable.  The 

results  of  using  this  finite  difference  scheme  to  solve  equations  (3.37)  and  (3.43)  will 
be  presented  in  Section  4.3.4. 

3.4-5  Steady  state  equation.  This  subsection  describes  the  steady  state  two 
dimensional  advection-diffusion  equation  and  the  numerical  method  used  to  solve  it 
in  this  research.  Following  reference  (25),  if  Ct  — +  0  as  t  — >  oo  and  kx  is  negligible  in 
equation  (3.37),  then 


uc. 


-  AryCyy  -j" 


(3.49) 
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with  initial  conditions 


c{0,y,z)  =  sin(ir(y  +  10)/20)  ■  8in(7r(2;  +  10)/20) 
c{x,y,z)=  0  y,z-^  ±10 


where  «  is  the  wind  speed  (positive  constant),  ky  =  both  equal  u  in  this  research 
and  are  the  y-aocis  and  z-axis  difFusivity  terms  respectively,  and  q  is  the  source  term. 
The  central-space  finite  difference  scheme  used  to  solve  equation  (3.49)  is 


u 


"Tn+l.p 


—  c 


m— l,p  _ 


2Ax 


(^v/Ay^)(C+i.p  -  2C,p  +  C-i.p) 

+  -  2c"  „  -I-  c"  p_i). 


(3.51) 


is 


The  solution  to  equation  (3.51)  will  be  compeired  to  the  exact  solution  which 


7r%  .  7r(y-f  10)^  .  ^n-(2-f- 10)^ 


20 


20 


(3.52) 


The  stability  condition  analysis  is  similar  to  the  two-dimensional  advection- 
diffusion  equation  discussed  in  Section  3.4.3. 


3.5  Conclusion 

Three  Gaussian  plume  models  are  described  in  this  chapter,  SCREEN,  AFTOX, 
and  GAUSPLUM.  SCREEN  uses  a  form  that  includes  reflection  terms,  but  that 
doesn’t  include  crosswind  terms.  SCREEN,  therefore,  mainly  does  calculations  of 
pollutant  concentrations  under  the  centerline  of  the  plume.  AFTOX  uses  three  forms 
of  the  Gaussian  model  which  include  two  Gaussian  puff  models,  one  when  there  is 
no  inversion  in  the  atmosphere,  and  one  when  there  is  an  inversion.  AFTOX  also 
uses  the  Gaussian  plume  model.  GAUSPLUM  uses  the  Gaussian  plume  model.  The 


results  and  comparison  of  SCREEN,  AFTOX,  the  progrcim  GAUSPLUM  developed 
in  this  study,  will  be  presented  in  Section  4.2.  The  other  models  of  air  pollution 
transport  discussed  in  this  chapter  all  follow  from  the  three-dimensional  ad vect ion- 
diffusion  equation.  They  are  all  special  cases  of  the  three-dimensional  advection- 
diffusion  equation.  The  advection  equation  does  not  include  any  diffusion.  The  one- 
«md  two-dimensional  advection-diffusion  equations  only  include  the  downwind,  and 
downwind  and  crosswind  diffusion  terms,  respectively.  The  steady-state  equation 
has  no  change  in  the  concentration  with  respect  to  time.  The  results  of  the  numer¬ 
ical  schemes  used  to  solve  the  advection  equation,  one,  two  and  three  dimensional 
advection-diffusion  equations,  and  the  two-  dimensional  steady  state  equation  will 
be  discussed  in  Sections  4.3.1,  4.3.2,  4.3.3,  4.3.4,  and  4.3.5  respectively. 
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IV.  Results 


4-1  Introduction 

The  purpose  of  this  study  is  to  examine  models  which  use  analytical  solutions 
in  their  calculations  of  pollutant  concentrations  and  to  find  numerical  solutions  of 
the  various  partial  differential  equations  used  in  pollution  modeling.  This  chapter 
includes  a  comparison  of  the  models  SCREEN  (10),  AFTOX  (19),  and  the  pro¬ 
gram  GAUSPLUM  (developed  in  this  study)  using  the  examples  taken  from  Turners 
workbook  (27)  described  in  Section  3.3.  The  results  of  the  numerical  solutions  of 
the  different  equations  described  in  Section  3.4  are  also  presented  here. 

The  comparison  of  SCREEN,  AFTOX,  and  GAUSPLUM  will  be  discussed  in 
Section  4.2  The  results  of  the  numerical  schemes  used  to  solve  the  partial  differential 
equations  and  the  sections  they  are  discussed  as  follows:  the  advection  equation  in 
Section  4.3.1,  the  one,  two  and  three  dimensional  advection-diffusion  equations  will 
be  discussed  in  Sections  4.3.2,  4.3.3,  and  4.3.4,  respectively,  and  two-dimensional 
steady  state  equation  in  Section  4.3.5. 

4.2  Comparison  of  Gaussian  models 

A  comparison  of  the  three  Gaussian  models  which  use  analytical  solutions 
in  their  calculations  is  done  using  the  six  problems  from  Turner’s  workbook  (27) 
described  in  Section  3.3.  Table  4.1  shows  the  comparison  of  SCREEN,  AFTOX, 
GAUSPLUM  using  problems  1,  2,  and  3,  which  calculate  the  concentration  at  a 
certain  location  where  a  receptor  is  located.  Table  4.2  shows  the  comparison  of 
SCREEN,  AFTOX,  GAUSPLUM  using  problems  4,5,  and  11,  which  give  the  location 
and  value  of  the  maximum  concentration.  The  models  all  have  solutions  to  the  six 
problems  within  an  order  of  magnitude  of  the  solutions  given  in  Turner’s  workbook. 
This  is  expected  since  they  each  use  a  form  of  the  Gaussian  plume  model. 


4-1 


Turner  Workbook  Examples  (27) 

ex. 

SCREEN 

AFTOX 

GAUSPLUM 

# 

1 

11.0  E-6 

11.34  E-6 

4.0  E-6 

10.95  E-6 

2 

33.0  E-6 

41.45  E-6 

23.0  E-6 

32.9  E-6 

3 

13.0  E-6 

NA 

9.0  E-6 

12.95  E-6 

Table  4.1  Comparison  of  concentrations  at  given  point 


As  Table  4.1  shows,  the  SCREEN  model  produces  higher  estimates  than  Turn¬ 
er’s  workbook  results  for  problems  1  and  2  and  isn’t  applicable  for  problem  3. 
This  result  is  because  the  equation  (see  eq.  3.1)  used  in  the  model  lacks  the  term 
exp[— found  in  both  AFTOX  and  GAUSSPLUM  which  allows  the  receptor 
to  be  at  points  away  from  the  x-ajcis.  The  AFTOX  model  consistently  underesti¬ 
mated  the  values  for  each  of  the  three  problems  in  Table  4.1.  All  three  models  results 
are  within  an  order  of  magnitude  on  either  side  of  the  results  in  Turners  workbook. 


The  results  are  similar  for  the  problems  which  determine  the  maximum  con¬ 
centration  and  its  location.  Table  4.2  shows  these  results.  Again  all  three  models 
solutions  are  within  an  order  of  magnitude  of  the  problems  in  Turners  workbook. 


Turner  Workbook  Examples  (27) 

ex. 

# 

Turner 
(gf/m^)  at  m 

SCREEN 
(g/m^)  at  m 

AFTOX 
{gjm^)  at  m 

GAUSPLUM 
{glm^)  at  m 

4 

280.0  E-6  at  1000 

235.6  E-6  at  1005 

180.0  E-6  at  1420 

263.4  E-6  at  1000 

5 

1.1  E-4  at  5600 

0.782  E-4  at  5454 

1.0  E-4  at  2839 

1.13  E-4  at  5475 

11 

6.4  E-5  at  13000 

2.52  E-5  at  12434 

6.8  E-5  at  6109 

6.08  E-5  at  12500 

Table  4.2  Comparison  of  maximum  concentration  calculations 


As  described  in  Chapter  III  both  the  SCREEN  (Section  3.2.1)  and  AFTOX 
(Section  3.2.2)  models  require  additional  input  such  as  ambient  temperature  (outside 
air  temperature)  and  exit  velocity  of  the  pollutant.  Figure  4.1  is  output  from  the 
SCREEN  program  and  Figure  4.2  is  from  the  AFTOX  program  for  problem  2  in 
Turners  workbook.  The  following  is  a  brief  description  of  the  two  figures.  Each 
program  asks  for  ambient  air  temperature,  emmission  rate,  wind  speed,  stack  height, 
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stack  gas  temperature,  and  the  location  of  the  receptor  for  the  calculation.  Problem 
2  says  it  is  an  overcast  winter  morning.  SCREEN  asks  for  the  ambient  temperature 
only,  while  AFTOX  asks  for  both  the  ambient  temperature  and  the  time  and  date 
of  the  calculation.  The  date  chosen  for  AFTOX  for  this  problem  is  lanuary  6,  1992. 
For  these  programs  the  value  of  the  ambient  temperature  is  taken  from  climatological 
data  for  the  WPAFB  area  for  the  time  of  year  described  in  each  problem.  The  exit 
velocity  and  stack  diameter  used  here  come  from  other  examples  in  the  workbook. 

The  program  GAUSPLUM  simply  asks  for  the  stability  cleiss,  emission  rate, 
wind  speed,  stack  height  and  the  location  of  the  receptor.  The  output  for  GAUS¬ 
PLUM  displays  the  value  of  the  distance  and  the  concentration  value.  It  only  does 
calculations  for  one  distance  at  a  time  so  for  problems  4,  5,  and  11  this  program  is 
run  several  times  to  determine  the  majcimum  concentration  and  its  distance. 
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SCREEN-1.1  MODEL  RUN 


♦**  VERSION  DATED  88300 

SIMPLE  TERRAIN  INPUTS: 

*** 

SOURCE  TYPE 

= 

POINT 

EMISSION  RATE  (G/S) 

= 

80.00 

STACK  HEIGHT  (M) 

= 

60.00 

STK  INSIDE  DIAM  (M) 

= 

2.00 

STK  EXIT  VELOCITY  (M/S) 

= 

2.00 

STK  GAS  EXIT  TEMP  (K) 

= 

293.50 

AMBIENT  AIR  TEMP  (K) 

= 

284.00 

RECEPTOR  HEIGHT  (M) 

= 

.00 

lOPT  (1=URB,2=RUR) 

= 

2 

BUILDING  HEIGHT  (M) 

= 

.00 

MIN  HORIZ  BLDG  DIM  (M) 

.00 

MAX  HORIZ  BLDG  DIM  (M) 

= 

.00 

BUOY.  FLUX  =  .03  H**4/S**3;  MOM.  FLUX  =  3.99  M**4/S**2. 

♦♦♦  STABILITY  CLASS  4  ONLY  *♦* 

***  10-METER  WIND  SPEED  OF  6.0  M/S  ONLY  *** 

***  SCREEN  DISCRETE  DISTANCES  *** 


CALCULATION  MAX  CONC  DIST  TO  TERRAIN 

PROCEDURE  (UG/M**3)  MAX  (M)  HT  (M) 


SIMPLE  TERRAIN  41.45  500.  0. 


Figure  4.1  Partial  output  of  SCREEN  program  for  problem  2 


^.3  Results  and  comparison  of  exact  solutions 

This  section  shows  the  results  of  the  finite  difference  methods  used  to  solve  the 
partial  differential  equations  described  in  chapter  III. 
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USAF  TOXIC  CHEMICAL  DISPERSION  MODEL 
AFTOX 

HPAFB  OH 
DATE:  01-06-1992 
TIME:  0800  LST 

CONTINUOUS  BUOYANT  PLUME 

TEMPERATURE  =  6  C 

HIND  DIRECTION  =  360 

HIND  SPEED  =  6  M/S 

SUN  ELEVATION  ANGLE  IS  5  DEGREES 

CLOUD  COVER  IS  8  EIGHTHS 

CLOUD  TYPE  IS  LOH  (St,  Ns,  FOG) 

GROUND  IS  DRY 

THERE  IS  NO  INVERSION 

ATMOSPHERIC  STABILITY  PARAMETER  IS  3.5 

EMISSION  RATE(KG/MIN)  »  4.8 

EFFLUENT  IS  STILL  BEING  EMITTED 

STACK  HEIGHT  ABOVE  GROUND (M)  =  60 

GAS  STACK  TEMP(C)  *  15 

VOLUME  FLOH  RATE(M3/MIN)  =  4 

EFFECTIVE  PLUME  HEIGHT (M)  =  60  AT  DISTANCE (M)  =  2 

CONCENTRATION  AVERAGING  TIME  IS  10  MIN 

HEIGHT  ABOVE  GROUND  IS  0  M 

DOHNHIND  DISTANCE  IS  500  M 

CROSSHIND  DISTANCE  IS  0  M 


THE  CONCENTRATION  IS  .018  MG  M-3 


Figure  4.2  Partial  output  of  AFTOX  program  for  problem  2 
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4-S.l  Advection  equation.  The  numerical  solution  to  the  advection  equa¬ 
tion  (see  eq.  3.11)  is  found  using  the  Lax-Friedrichs  and  the  leapfrog  finite  difference 
schemes.  For  both  schemes  A  =  0.8,  Ax  =  0.1,  and  At  =  0.08.  On  the  boundary, 
when  X  =  —2,  c  =  0,  and  when  m  =  3.  The  finite  difference  for¬ 

mula,  equation  (3.17),  in  Section  3.4.1  is  used  in  the  Ada  program  in  appendix  (B.l) 
in  procedure  ComputeJNew.  Figure  4.3  shows  the  solution  to  the  initial-boundary 
value  problem  (see  eq.  3.12)  discussed  in  Section  3.4.1  using  the  Lax-Friedrichs  finite 
difference  scheme.  Figure  4.4  shows  the  solution  to  the  same  initial- boundary  value 
problem  using  the  leapfrog  finite  difference  scheme,  equation  (3.18),  in  procedure 
Compute_New. 


X 


Figure  4.3  Lax-Friedrichs  solution  of  advection  equation 

In  both  Figure  4.3  and  Figure  4.4  the  exact  solution  is  shown  as  a  dashed  line, 
and  the  solution  of  the  finite  difference  scheme  is  shown  as  the  curve  with  diamonds. 
The  leapfrog  method  has  more  oscillations  in  its  solution  than  the  Lax-Friedrichs, 
but  the  overall  accuracy  is  better  for  the  leapfrog  scheme.  The  approximation  at 
the  peak  in  Figure  4.4,  the  leapfrog  solution,  is  much  better  than  in  Figure  4.3.  The 
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Figure  4.4  Leapfrog  solution  of  advection  equation 

solution  of  the  Lax-Friedrichs  scheme  can  be  improved  by  decreasing  the  value  of 
Ax  while  keeping  the  same  value  of  A. 

As  discussed  in  Section  3.4.1  the  Lax-Friedrichs  scheme  is  conditionally  stable. 
To  show  this  A  =  1.6  is  used  in  the  scheme  with  the  results  shown  in  Figure  4.5. 

In  Figure  4.5  the  exact  solution  is  shown  as  a  solid  line  while  the  solution  of  the 
Lax-Friedrichs  scheme  is  shown  as  a  line  with  diamonds. 

Figures  4.3  and  4.4  show  that  the  solutions  of  the  Lax-Friedrichs  scheme  and 
the  leapfrog  scheme  are  reasonable  approximations  to  the  solution  of  the  advection 
equation.  As  the  values  of  Ax  and  At  are  decreased,  while  keeping  A  constant,  the 
solutions  of  the  schemes  become  better  approximations  to  the  advection  equation. 
In  the  next  section  a  diffusion  term  is  added  to  the  advection  equation  and  the 
one-dimensional  advection-diffusion  equation  is  looked  at. 
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Figure  4.5  Laix- Friedrichs  solution  with  A  =  1.6 

4.S.S  1-D  advection-diffusion  equation.  The  one-dimensional  advection- 
diffusion  equation  (see  eq.  3.22)  is  a  combination  of  the  advection  equation  described 
in  Section  4.3.1  and  the  diffusion  equation  discussed  in  Appendix  C.l.  The  forward¬ 
time  central-space  scheme  (see  eq.  3.23)  is  used  to  solve  the  advection-diffusion 
equation  (see  eq.  3.22).  It  is  in  Ada  procedure  Compute_New  which  is  found  in 
Appendix  B.2. 

The  stability  condition  for  this  method,  discussed  in  Section  3.4.2,  is  <1/2 
where  mu  =  AtfAx^.  This  stability  condition  means  the  time  step  At  is  at  most 
Ax^f2kx  and  when  decreasing  Ax  by  half  to  increase  the  spatial  accuracy,  At  must 
decrease  by  one-fourth.  This  restriction  limits  practical  use  of  this  method.  The 
data  in  Table  4.3  show  the  cases  used  in  this  method. 

Table  4.3  includes;  u,  the  advection  coefficient;  kx,  the  diffusion  coefficient; 
Ax,  the  space  incrementer;  At,  the  time  incrementer;  the  range  in  space;  kxfi,  which 
from  the  stability  analysis  must  be  less  than  or  equal  to  one  half;  a  =  {Axu)/{2kx), 
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Try 

# 

u 

■ 

Ax 

At 

X- range 

a  =  (Axu)J{2kx) 
(<  1-0) 

t 

1 

4 

■a 

1 

0.1 

-10.. 10 

0.4 

0.5 

2.0 

2 

10 

4 

1 

0.1 

-10.. 10 

0.4 

1.25 

2.0 

3 

8.1 

4 

1 

0.1 

-10..10 

0.4 

1.0125 

2.0 

4 

6 

7 

2 

0.2 

-20..20 

0.35 

0.85714 

10.0 

5 

1 

1 

0.5 

0.05 

-10..10 

0.2 

0.25 

2.0 

6 

3.95 

1 

-10..10 

0.2 

0.9875 

2.0 

mm 

4.01 

1 

0.5 

0.05 

-10..10 

0.2 

1.0025 

2.0 

8 

1.0025 

0.25 

0.5 

0.05 

-10..10 

0.05 

1.0025 

2.0 

9 

1.0025 

0.25 

0.5 

0.2 

-10..10 

0.2 

1.0025 

2.0 

10 

4 

0.25 

0.5 

0.25 

-10..10 

0.15 

4 

2.0 

Table  4.3  1-D  advection-diffusion  test  data 


which  Strikwerda  shows  must  be  less  than  or  equal  to  one;  and  t,  the  time  of  final 
calculation. 

The  cases  included  in  Table  4.3  all  satisfy  the  stability  condition  kxfi  <  1/2. 
These  cases  were  done  to  test  the  other  condition  a  =  {Axu)f{2kx)  <  1  to  see  if 
oscillations  detracted  from  the  accuracy  of  the  solution.  The  cases  also  varied  the 
advection  term  to  see  its  impact  on  the  solutions.  In  Figure  4.6  the  solution  for  test 
case  number  6  is  plotted. 

In  Figure  4.7  the  solution  for  test  case  number  10  is  plotted.  Case  number 
10  has  an  oscillating  solution  because  the  a  <  1  condition  is  not  met.  In  this  case 
a  =  4.0  which  causes  the  oscillation  even  though  the  stability  condition  is  met. 

These  two  cases,  numbers  6  and  10  in  Table  4.3,  have  similar  advection  terms, 
3.95  and  4.0,  respectively.  The  diffusion  terms  differ  by  a  larer  margin  and  contribute 
to  the  differences  in  the  a  condition  and  also  to  the  oscillations  in  case  10.  Case 
6  satisfies  the  a-condition  and  has  a  non-oscillating  solution,  while  case  10  violates 
the  condition  and  has  an  oscillating  solution. 

In  the  other  cases  of  Table  4.3  the  advection  coefficient  and  the  diffusion  coef¬ 
ficient  are  varied  to  see  the  effect  of  a  higher  advection  term  than  diffusion  term.  In 
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all  of  the  cases  with  a  greater  than  one  the  solution  oscillates  as  expected.  In  the 
other  cases,  where  q  is  less  than  one  and  the  advection  coefficient  is  greater  than 
the  diffusion  coefficient,  the  solution  is  similar  to  case  number  six  in  Figure  4.6.  In 
the  next  section,  a  second  diffusion  coefficient  is  added  to  the  equation  to  see  if  the 
diffusion  has  a  larger  impact  on  the  solution. 


4. 3. 3  2-D  advection-diffusion  equation.  The  two-dimensional  advection- 
diffusion  equation  (see  eq.  3.29)  adds  the  crosswind  (y-axis)  diffusion  term  to  the 
one- dimensional  equation.  The  forward-time  central-space  scheme,  equation  (3.30), 
in  Section  3.4.3  is  used  to  solve  the  2-D  advection-diffusion  equation  (3.29).  It  is 
used  in  Ada  procedure  Compute-Final  found  in  Appendix  B.3.  The  data  in  Table  4.4 
show  the  cases  used  in  this  method. 


H 

u 

ICx 

HI 

Ax 

Ay 

At 

time 

1 

2 

0.0625 

0.0625 

0.2 

0.2 

0.02 

1.6 

2 

0.5 

0.0625 

0.0625 

0.2 

0.2 

0.02 

1.6 

3 

4 

0.05 

0.05 

0.2 

0.2 

0.02 

1.6 

4 

a 

0.05 

0.05 

0.1 

0.1 

0.08 

1.6 

5 

0.5 

0.5 

0.5 

0.2 

0.2 

0.02 

1.6 

6 

0.5 

1 

1 

0.2 

0.2 

0.1 

1.6 

mm 

0 

1 

1 

0.2 

0.2 

0.1 

1.6 

8 

0 

1 

1 

0.2 

0.2 

0.08 

1.6 

9 

0.5 

1 

1 

0.2 

0.2 

0.04 

1.6 

10 

0.5 

0.9 

0.9 

0.2 

0.2 

0.04 

1.6 

11 

0.5 

1 

1 

0.2 

0.2 

0.01 

1.6 

12 

0.5 

1 

1 

0.2 

0.2 

0.02 

1.6 

13 

1 

1 

1 

0.2 

0.2 

0.01 

1.6 

14 

4 

1 

1 

0.2 

0.2 

0.01 

1.6 

15 

0 

1 

1 

0.2 

0.2 

0.01 

1.6 

16 

0 

0.0625 

0.0625 

0.2 

0.2 

0.02 

1.6 

17 

1 

0.0625 

0.0625 

0.2 

0.2 

0.02 

1.6 

18 

1.5 

0.0625 

0.0625 

0.2 

0.2 

0.02 

1.6 

Table  4.4  2-D  advection-diffusion  test  data 
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Table  4.4  includes;  u,  the  advection  coefficient;  k^,  ky,  the  x-axis,  and  y-axis 
diffusion  coefficients;  Ax,  the  x*axis  incrementer;  Ay,  the  y-axis  lacrementer;  At, 
the  time  incrementer;  and  t,  the  time  of  final  calculation. 

The  stability  condition  for  this  method,  discussed  in  Section  3.4.3,  is  kj.fi  <  1/4 
where  ft  =  At/Ax^.  This  stability  condition  means  the  time  step  At  is  at  most 
Ax^f2kx  and  when  decreasing  Ax  by  half  to  increase  the  spatial  accuracy.  At  must 
decrease  by  one-eighth.  This  restriction  limits  practical  use  of  this  method  even 
more  than  the  one-dimensional  case. 

Two  cases  that  show  the  difference  between  complying  with  the  stability  condi¬ 
tion  and  violating  it  are  case  5  and  case  12.  Case  5  is  an  example  of  complying  with 
the  stability  condition,  while  case  12  is  one  that  violates  the  condition.  The  only 
differences  between  the  two  cases  are  that  for  case  5,  kj  =  1/2  and  kjAtlA^x^  =  1/4, 
while  for  case  12,  fci  =  1.0  and  kxAtfAx^  =  1/2. 

Figure  4.8  shows  the  solution  to  equation  (3.32)  using  test  case  5  from  Table  4.4, 
where  the  view  is  looking  at  the  xz- plane  with  the  x-axis  going  from  left  to  right, 
and  the  concentration  is  on  the  z-axis.  The  y-axis  goes  into  the  page.  The  lines 
in  the  figure  are  the  concentration  contours  for  different  values  of  y  similar  to  that 
in  Figure  4.9  which  shows  the  solution  using  case  5  v/ith  y  =  0.  When  y  =  0  the 
resulting  contour  is  the  maximum  concentration  contour  since  it  is  directly  downwind 
from  the  source.  The  minimum  concentration  contour  is  near  the  boundary  at  y  = 
±10.  Figure  4.10  shows  the  solution  to  equation  (3.32)  using  test  case  12  with  the 
same  view  as  Figure  4.8.  Case  12  violates  the  stability  condition  which  results  in 
large  oscillations.  Again  the  lines  in  the  figure  are  the  concentration  contours  for  a 
given  value  of  y.  Figure  4.11  shows  the  solution  using  case  12  with  y  =  0  and  shows 
the  maximum  concentration  contour. 

Most  of  the  cases  in  Table  4.4  vary  the  terms  u,  kx,  ky,  and  time  while  keeping 
the  Ax  and  Ay  constant.  One  exception  to  this  is  case  3  and  case  4  which  changes 
Ax,  Ay  and  At,  but  keeps  the  other  terms  constant.  Case  3  is  stable  since  kjfi  =  1/4 


4-12 


which  meets  the  stability  condition.  Case  4  on  the  other  hand  does  not  meet  the 
stability  condition  and  is,  therefore,  unstable. 


c(x.y) 


-2  -1.5  -1  -0.5  0  0.5  1  1.5  2  2.5 


X 

Figure  4.8  2-D  test  case  5 

These  results  show  how  only  slightly  violating  the  stability  condition  can  cause 
large  oscillations.  The  same  result  will  be  shown  in  the  next  section  for  the  3-D 
advection-diffusion  equation. 


3 


Figure  4.11  2-D  test  case  12  with  y  =  0 

4.3.4  3~D  advection- diffusion  equation.  The  three-dimensional  advec- 

tion-diffusion  equation  (3.37)  adds  the  vertical  (z-axis)  diffusion  term  to  the  two- 
dimensional  equation.  The  other  equations  are  all  special  cases  of  the  three-dimen¬ 
sional  advection-diffusion  equation.  The  advection  equation  has  no  diffusion  term. 
The  one-dimensional  advection-diffusion  equation  has  no  crosswind  or  vertical  diffu¬ 
sion,  and  the  two-dimensional  advection-diffusion  equation  has  no  vertical  diffusion. 
The  steady-state  equation  has  no  concentration  change  with  respect  to  time.  This 
research  looks  at  the  three-dimensional  advection-diffusion  equation  in  two  ways, 
one  with  the  source  as  part  of  the  initial  conditions,  and  a  second  way  by  adding  a 
forcing  function  or  source  term  so  that  initial  conditions  are  incorporated  into  the 
equation.  This  section  will  look  at  the  results  of  using  a  forward-time  central-space 
scheme,  equation  (3.40)  in  Section  3.4.4,  to  solve  the  three-dimensional  advection- 
diffusion  equation  with  the  forcing  function  added,  equation  (3.43).  The  scheme  is 
used  in  Ada  procedure  Compute  Jinal  found  in  Appendix  B.4.  Table  4.5  shows  the 
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cases  used  in  the  method  without  the  source  term  and  with  the  initial  conditions, 
equation  (3.38),  discussed  in  Section  3.4.4. 


Case 

# 

u 

kx 

At 

A* 

stable 
<  0.125 

t 

1 

0.5 

1 

1 

2 

0.25 

10 

2 

0.5 

0.5 

1 

2 

0.125 

10 

3 

0.5 

0.25 

1 

2 

0.0625 

10 

4 

0.5 

0.125 

1 

2 

0.03125 

10 

Table  4.5  3-D  without  source  term 


In  Table  4.5;  u  is  the  advection  coefficient,  kx  =  ky  =  k^  are  the  diffusivity  coef¬ 
ficients,  At  is  the  change  in  time,  Ai  is  the  change  in  the  spatial  directions  where 
i  is  X,  y,  or  z,  stable  is  the  stability  conditon,  and  t  is  the  time  of  the  calculation. 
In  these  cases  the  only  difference  between  them  is  the  change  in  At  which  affects 
the  stability  condition.  The  three  previous  sections  show  that  the  o-condition  must 
be  met  for  the  solution  to  be  well-behaved.  In  this  section  the  stability  condition  is 
examined,  that  is  kxAtjAi^  <  0.125,  where  i  is  x,  y,  or  z. 

Figure  4.12  shows  case  3  from  Table  4.5,  as  a  line  with  squares,  compared  to 
the  exact  solution,  equation  (3.42),  as  a  solid  line.  The  figure  shows  the  solution  for 
X  =  8,  2  =  0,  and  t  =  10. 

Figure  4.13  shows  cases  1,  2  and  3  from  Table  4.5,  as  a  line  with  crosses,  a 
line  with  diamonds  and  a  line  with  squares,  respectively,  compared  to  the  exact 
solution,  equation  (3.42),  as  a  solid  line.  The  figure  again  shows  the  solution  for 
X  =  8,  2  =  0,  and  i  =  10. 

The  reason  the  numerical  solution  only  seems  to  move  toward  the  exact  solution 
near  y  =  0  and  not  change  much  near  the  boundary  is  that  the  boundary  condition 
c(x,i/,z,t)  =  0  when  x  =  ±10  or  y  =  ±10,  or  2  =  ±10  does  not  represent  what 
the  exact  solution  does  at  those  boundary  points.  The  exact  solution  has  the  same 
boundary  condition,  except  that  the  boundary  is  at  ±00.  The  boundary  condition 
for  the  numerical  solution  causes  the  boundary  to  act  as  a  sink  for  the  pollutant 
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c{8.y.0.10)  c{8.y.0.10) 
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as  it  is  dispersed  ajid  advected.  The  other  variation  described  in  Section  3.4.4  has 
conditions  such  that  the  exact  calculation  has  the  same  boundary  conditions  as  the 
numerical  calculation. 

Table  4.6  shows  the  cases  used  in  the  scheme  with  the  source  term,  equation 
(3.45). 


Case 

# 

u 

— 

kx 

Ai 

stable 
<  0.125 

At 

t 

1 

0.2 

1 

5 

0.04 

1 

10 

2 

0.2 

1 

5 

0.02 

0.5 

10 

3 

0.2 

1 

5 

0.01 

0.25 

10 

4 

0.2 

1 

5 

0.005 

0.125 

10 

Table  4.6  3-D  with  source  term 

In  Table  4.6:  u  is  the  advection  coefficient,  kx  =  ky  =  are  the  diffusivity  coeffi¬ 
cients,  Ai  is  the  change  in  the  spatial  directions  where  i  is  x,  y,  or  z,  stable  is  the 
stability  conditon.  At  is  the  change  in  time,  and  t  is  the  time  of  the  calculation. 
Again,  At  is  the  only  difference  between  the  cases  and  this  changes  the  stability 
condition  as  shown  in  the  table.  These  cases  all  satisfy  the  stability  condition.  The 
following  three  figures  are  of  cases  1,  2,  and  3  which  show  that  <is  At  is  decreased  the 
numerical  solution  converges  on  the  exact  solution.  This  result  would  be  the  same 
as  decreasing  both  the  spatial  terms  Ai  and  At  such  that  the  stability  condition 
remained  the  same  as  in  Table  4.6.  Figure  4.14  shows  case  1  as  a  line  with  squares. 
Figure  4.15  shows  case  2  as  a  line  with  crosses.  Figure  4.16  shows  case  3  as  a  line 
with  diamonds.  In  all  three  figures  the  exact  solution,  equation  (3.47),  is  a  solid  line. 

The  tvvo  variations  of  the  three-dimensional  advection-diffusion  equation  in 
this  section  show  how  the  boundary  conditions  can  have  an  effect  on  the  numerical 
solution.  The  boundary  condition  for  the  variation  without  the  source  term  does 
not  represent  the  actual  conditions  where  the  boundary  is  at  ±cx).  The  numerical 
solution  is  less  than  actual  since  the  boundary  condition  artificially  lowered  the  value 
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c(5.y,0.10) 


Figure  4.16  3-D  case  3  with  source  term 

at  points  near  the  boundary.  The  other  variation  has  the  same  boundary  conditions 
as  the  exact  calculation  and  the  numerical  and  exact  solutions  are  much  closer. 

4.3.5  Steady  state  equation.  The  two-dimensional  steady  state  equation  is 
similar  to  the  two-dimensional  advection-difFusion  equation  except  the  concentration 
does  not  change  with  respect  to  time.  This  section  will  look  at  the  results  of  using  a 
central-space  scheme,  equation  (3.51)  in  Section  3.4.5,  to  solve  the  two-dimensional 
steady  state  equation  (3.49)  with  the  forcing  function  added  in  the  initial  conditions, 
equation  (3.50).  The  scheme  is  used  in  Ada  procedure  Compute_Final  found  in 
Appendix  B.5.  Table  4.7  shows  the  cases  used  in  this  scheme. 

In  Table  4.7  u  is  the  advection  coefficient,  ky  =  =  u  are  the  diffusivity  coefficients, 

Ax  is  the  change  in  x-axis  direction,  Ai  is  the  change  in  the  crosswind  and  vertical 
directions  where  i  is  y,  or  z,  and  stable  is  the  stability  conditon  {kyAx)/{Ay^).  In 
these  cases  the  only  difference  between  them  is  the  change  in  Ax  which  also  changes 
value  of  the  stability  condition  by  the  same  factor. 
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Case 

# 

1 

Ax 

Ai 

stable 
<  0.25 

1 

a 

n 

5 

5 

0.2 

2 

D 

D 

2 

5 

0.08 

3 

a 

11 

1 

5 

0.04 

4 

a 

a 

0.5 

5 

0.02 

Table  4.7  2-D  steady  state  data 


Figure  4.17  shows  case  1  as  a  line  with  squares.  Figure  4.18  shows  case  3  as 
a  line  with  crosses.  Figure  4.19  shows  case  4  as  a  line  with  diamonds.  Figure  4.20 
shows  a  comparison  of  all  four  cases  with  the  exact  solution.  In  all  four  figures  the 
exact  solution,  equation  (3.52),  is  a  solid  line. 


Figure  4.17  Steady  state  case  1 


If  the  stability  condition  is  satisfied,  then  decreasing  the  value  of  Ax  does  not 
require  any  changes  in  other  variables.  However,  if  the  stability  condition  is  satisfied, 
and  a  decrease  in  Ay  and  Az  is  wanted  by  one  half,  then  Az  must  be  decreased  by 
one  fourth. 
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c(10.y,0)  c(10.y.0) 


y 

Figure  4.18  Steady  state  case  3 
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Figure  4.20  Steady  state  comparison 

4-4  Conclusion 

The  results  of  the  comparison  of  SCREEN,  AFTOX,  eind  GAUSPLUM  show 
that  they  can  be  used  in  general  cases  such  as  those  in  the  problems  of  Turner’s 
workbook,  SCREEN  and  AFTOX  can  also  be  used  in  much  more  detailed  ceises 
since  they  each  ask  for  more  details  than  provided  by  the  problems  in  the  workbook. 

The  numerical  solutions  found  for  the  advection  equation,  all  three  advection- 
diffusion  equations,  and  the  steady-state  equation  show  the  effect  of  violating  the 
stability  condition  or  the  condition  that  handles  oscillations  in  the  solution. 

Each  of  these  results  show  possibilities  for  further  research  in  the  area  of  air 
pollution  modeling,  some  of  which  will  be  discussed  in  Chapter  V. 
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V.  Conclusion  and  Recommendations 

5. 1  Conclusion 

The  purpose  of  this  study  was  to  analyze  models  which  use  analytical  solutions 
in  their  calculations  of  pollutant  concentrations  and  to  find  numerical  solutions  of 
the  various  partial  differential  equations  used  in  pollution  modeling. 

Many  air  pollution  transport  models  use  Gaussian  analytical  equations  to  solve 
the  partial  differential  equations  used  in  pollution  modeling.  This  study  has  looked 
at  three  such  models  and  made  comparisons  between  them  using  example  problems 
which  deal  with  finding  pollutant  concentrations  at  a  given  location  and  time  and 
finding  the  maximum  concentration  and  its  location. 

The  Gaussian  plume  model  used  in  many  air  pollution  transport  models  has 
difficulties  when  the  wind  speeds  are  low  to  calm.  This  research  has  looked  at 
numerically  solving  the  advection  equation,  advection-diffusion  equation,  and  the 
steady-state  equation  using  low  wind  speeds  in  the  calculations.  The  results  of  these 
calculations  were  compared  to  exact  solutions. 

5.2  Recommendations 

Pollution  modeling  has  many  areas  of  interest  to  consider  for  further  research 
based  on  this  research.  The  first  is  to  look  at  more  of  the  models  which  use  analytical 
equations  in  their  calculations.  This  research  looked  at  SCREEN,  AFTOX,  and 
GAUSPLUM.  Two  other  models  to  consider  are  TOXST  (12),  and  TOXLT  (11) 
which  are  referenced  in  the  bibliography  of  this  thesis.  These  models  look  at  short 
term  and  long  term  exposure  of  pollutants. 

Another  area  of  possible  further  research  could  be  continuing  the  two  and  three- 
dimensional  research  by  looking  at  variable  advection  and  diffusion  coefficients.  This 
research  looked  at  those  coefficients  as  constant  in  determining  solutions.  Although 
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the  advection  coefficient  is  often  averaged  as  the  mean  wind  speed,  it  can  be  a  made 
variable  function  dependent  on  some  standard  deviation  of  the  wind  speed.  The 
diffusion  coefficient  is  in  reality  variable  and  dependent  on  the  distance  from  the 
source.  These  variations  could  be  added  to  numerical  (finite-  difference)  schemes. 
Also,  other  numerical  schemes,  such  as  the  finite  element  method,  could  be  applied 
to  the  same  equations  evaluated  in  this  research. 

Another  subject  of  possible  further  research  could  be  expanding  the  regions  of 
concern  and  using  parallel  computing  systems  in  finding  solutions  to  the  equations 
to  better  simulate  the  boundaries  in  real  life,  like  x,y,z  ±cx).  This  research  some¬ 
what  limited  in  range  due  to  storage  limits  encountered  using  the  Ada  progamming 
language. 

A  fourth  area  of  possible  further  research  could  be  looking  at  the  area  of  risk 
assessment  of  pollutants  done  by  models  currently  available  through  the  Environ¬ 
mental  Protection  Agency  and  other  sources  of  pollution  models. 
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Appendix  A.  GAUSPLUM  Ada  Code 


This  Appendix  contains  the  Ada  source  code  for  the  program  GAUSPLUM. 
The  description  for  the  variables  used  in  the  programs  is  within  each  program  in 
the  form  of  comment  blocks  before  each  procedure  within  the  program.  Comment 
lines  start  with  a  double  hyphen  and  comment  blocks  begin  and  end  with  a  line  of 
hyphens. 


—  FILE :  gausplum . a 

—  PROJECT:  Gaussian  Plume  Model 

—  DATE:  16  APR  93 

—  VERSION:  Version  1.0 

—  AUTHORS:  Capt  Dave  Paal 

—  DESCRIPTION:  This  program  calculates  the  concentration  of 

—  pollutants  in  a  plume  using  a  Gaussian  Plume  Equation.  The 

—  program  asks  the  user  for  the  source  strength,  the  surface 

—  wind  speed,  coordinates  of  the  receptor (x, y, z) ,  Pasquill 

—  Gifford  stability  condition,  and  the  effective  height  of  the 

—  plume.  The  diffusion  terms,  sigmas,  aore  calculated  using  the 

—  Pasquill  Gifford  sigmas,  and  x  given  by  the  user. 

—  OPERATING  SYSTEM:  UNIX/Sun  Sparc  Station 

—  LANGUAGE:  Meridian  Ada 
"  FILES  USED: 


—  CONTEXT  CLAUSES 

with  text_io; 
with  Math.Lib; 
use  Math_Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  gausmod  is 
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—  TYPE  DECLARATIONS 


—  type  Stable.Type  is  (A,B,C,D,E,F) ; 

—  GLOBAL  VARIABLES  AND  EXCEPTIONS 

q  :  float;  — Source  Strength,  grzuns/second 

X  :  float;  — Distance  in  X-direction  in  meters 

U  :  float;  — Surface  vind  speed  in  meters/second 

Y  :  float;  — Distance  in  Y-direction  in  meters 

Z  :  float;  — Distance  in  Z-direction  in  meters 

H  :  float;  — Effective  height  of  pl\ime  in  meters 

Stable  :  character;  — Stability  condition  A  to  F 
SigmaY  :  float;  — lateral  diffusion  term 

SigmaZ  :  float;  — vertical  diffusion  term 

Concentration  :  float;  — calculated  concentration 

—  PROCEDURE:  Get.Info 

—  DESCRIPTION:  This  procedure  asks  the  user  for  the  following: 

—  Source  strength,  coordinates  of  receptor  (X,Y,Z),  surface 

—  wind  speed,  stability  condition,  and  the  effective  height  of 

—  the  plume . 

—  INPUT  PARAMETERS:  Q  :  source  strenght  in  grams  per  second, 

X  :  distance  in  x-direction  in  meters, 

Y  :  distance  in  y-direction  in  meters, 

Z  :  distanct  in  z-direction  in  meters, 

U  :  surface  wind  speed  in  meters  per  second, 
H  :  effective  height  of  plume  in  meters. 
Stable  :  stability  condition,  character  A  to  F. 

—  OUTPUT  PARAMETERS:  Same  as  Input  Parameters. 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED:  Same  as  parameters. 

—  CALLED  BY:  main. 

—  CALLS :  None . 

procedure  Get_Info  (Q  :  in  out  float; 

X  :  in  out  float; 

Y  :  in  out  float; 

Z  :  in  out  float; 

U  :  in  out  float; 


A-2 


H  :  in  out  float; 

Stable  :  in  out  character)  is 

begin 

Text_Io .put_line("Enter  the  Source  Strength  in  grams  per 
second. ") ; 

My_Float_Io . get (Q) ; 

Text.Io .put_line("Enter  the  distance  in  the  X-direction  in 
meters . ; 

My_Float_Io.get(X) ; 

Text.Io .put_line ("Enter  the  Y-direction  distance  in 
meters.") ; 

My_Float_Io.get(Y) ; 

Text_Io .put_line("Enter  the  Z-direction  disteoice  in 
meters.") ; 

My_Float_Io.get(Z) ; 

Text _ I o. put .line ("Enter  the  surface  wind  speed  in 
meters/second.") ; 

My_Float_Io.get(U) ; 

Text.Io .put.line ("Enter  the  effective  height  of  the  plume  in 
meters.") ; 

My_Float.Io.get (H) ; 

Text.Io.put.i ineC'Enter  the  stability  condition,  in  upper 
case.") ; 

Text.Io. get (Stable) ; 
end  Get.Info; 


—  PROCEDURE:  Calc.PG.Sigma 

—  DESCRIPTION:  This  procedure  calculates  the  Pasquill  Gifford 

diffusion  terms  sigma-y,  auid  sigma-z  for  the  given  stability 
condition,  A  to  F,  and  the  disteuice  X  in  meters  from  the 
source . 

The  general  form  of  the  equations  are: 

sigmay  =  (K1  ♦  X)/[1.0  +  (X/K2)]**K3, 
and  sigmaz  =  (K4  ♦  X)/[1.0  +  (X/K2)]**K5. 

As  found  in  Air  Pollution  Modeling  by  Paolo  Zannetti  on 
p.l49.  The  constants  K1  to  K5  used  in  this  procedure  were 
derived  by  Gifford  in  1976  in  a  diffusion  experiment  in  flat 
terrain. 

—  INPUT  PARAMETERS:  Stable  :  stability  condition  given  by  user, 

X  :  x-direction  distance  of  receptor. 
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—  OUTPUT  PARAMETERS:  SigmaY  :  lateral  diffusion  term, 

SigmaZ  :  vertical  diffusion  term. 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED :  Same  as  parauneters . 

—  CALLED  BY:  Main. 

—  CALLS :  None . 


procedure  Calc_PG_Sigma  (Stable  :  in  character; 

X  :  in  Float; 

SigmaY  :  in  out  Float; 

SigmaZ  :  in  out  Float)  is 

begin 

case  Stable  is 
when  'A'  => 

SigmaY  :=  (X  ♦  0.250)/exp(0.189*ln(1.0  +  (X/927.0))); 
SigmaZ  :=  (X  *  0 . 1020)/exp(-l .918*ln(l .0  +  (X/927.0))); 
when  'B'  => 

SigmaY  :=  (X  ♦  0.202)/exp(0.162*ln(1.0  +  (X/370.0))); 
SigmaZ  :=  (X  ♦  0.0962)/exp(-0.101*ln(1.0  +  (X/370.0))); 
when  'C’=> 

SigmaY  :*  (X  *  0.134)/exp(0.134*ln(1.0  +  (X/283.0))); 
SigmaZ  :*  (X  *  0.0722)/exp(0.102*ln(1.0  +  (X/283.0))); 
when  ’D'  => 

SigmaY  :=  (X  ♦  0.0787)/exp(0.135*ln(1.0  +  (X/707.0))); 
SigmaZ  :=  (X  *  0.0475)/exp(0.465*ln(1.0  +  (X/707.0))); 
when  'E'  *> 

SigmaY  :=  (X  ♦  0 .0566)/exp(0. 137*ln(l .0  +  (X/1070.0))) ; 
SigmaZ  :=  (X  ♦  0.0335)/exp(0.624*ln(l.0  +  (X/1070.0))); 
when  'F'  => 

SigmaY  :=  (X  *  0.037) /exp (0.134* In (1.0  +  (X/1170.0))) ; 
SigmaZ  :=  (X  *  0 .022)/exp(0.7*ln(l .0  +  (X/1170.0))) ; 
when  others  => 

Text_Io .put_line("Invalid  stability  condition"); 
end  case; 

end  Cal c_PG_ Sigma; 

—  PROCEDURE:  Calc_Concentration 

—  DESCRIPTION:  This  procedure  calculates  the  pollution 

concentration  for  the  user  entered  source  strength, 
distance,  stability  condition  surface  wind  speed,  and 
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effective  height  of  the  plume. 

—  INPUT  PARAMETERS:  SigmaY  :  lateral  diffusion  term, 

SigmaZ  :  vertical  diffusion  term, 

Q  :  source  strength, 

H  :  effective  height  of  plume, 

Y  :  distance  in  y-direction, 

Z  :  distance  in  z-direction, 

U  :  surface  wind  speed. 

—  OUTPUT  PARAMETERS:  Concentration  :  calculated  concentration. 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED:  Same  as  parameters. 

—  CALLED  BY:  Main. 

—  CALLS :  None . 

procedure  Calc_Concentration  (SigmaY  :  in  Float; 

SigmaiZ  :  in  Float; 

Q  :  in  Float; 

Y  :  in  Float; 

Z  :  in  Float; 

U  :  in  Float; 

H  :  in  Float; 

Concentration  :  out  Float)  is 

result  1  :  float; 
resultZ  :  float; 
results  :  float ; 
result4  :  float; 

Pi  :  float  :=  3.1415926; 

begin 

result 1  :=  (Q/(2.0  *  Pi  ♦  SigmaY  ♦  SigmaZ  *  U)); 
result2  :=  exp((-l .0/2.0)*(Y/SigmaY)*(Y/SigmaY)) ; 
results  :=  exp((-l .0/2.0) ♦((Z-H) /SigmaZ) ♦((Z-H) /SigmaZ) ) ; 
result4  :=  esp((-1.0/2.0)*((Z+H)/SigmaZ)*((Z+H)/SigmziZ)) ; 
Concentration  :=  resultl  ♦  result2  *  (results  +  result4) ; 

end  Calc_Concentration; 


—  PROCEDURE:  Print_Concentration 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 
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It  prints  out  the  concentration  and  the  coordinates  of  the 
receptor,  the  source  strength,  effective  height  of  the  plume, 
surface  wind  speed,  and  the  stability  condition  given. 


—  INPUT  PARAMETERS: 

Concentration 

:  calculated  concentration. 

— 

X,Y,Z 

:  recptor  coordinates. 

— 

Q 

:  source  strength. 

— 

H 

:  effective  height  of  plume 

— 

U 

:  surface  wind  speed. 

— 

Stable 

;  stability  condition. 

—  OUTPUT  PARAMETERS: 

None. 

"  LOCAL  VARIABLES: 

None. 

—  GLOBALS  USED:  All 

• 

—  CALLED  BY:  Main. 

—  CALLS :  None . 

procedure  Print_Info  (Concentration 

X,  Y,  Z 

q 


u 

H 

stable 


in  Float; 
in  Float; 
in  Float; 
in  Float; 
in  Float; 
in  character)  is 


begin 

Text_Io.put_line("The  following  concentration  results 
from") ; 

Text.Io .put_line("these  input  variables:"); 
Text_Io.put_line("Source  strength  =  "); 

My_Float_Io.put(Q) ; 

Text_Io .new.line; 

Text_Io .put_line("X  distance  =  ") ; 

My_Float_Io . put (X) ; 

Text.Io . new.line ; 

Text.Io .put_line("Y  distance  =  ") ; 

My_Float_Io.put(Y) ; 

Text.Io .new.line; 

Text_Io.put_line("Z  distance  =  "); 

My_Float_Io . put (Z) ; 

Text.Io . new.line ; 

Text.Io .put_line("Surface  wind  direction  is  ") ; 
My_Float_Io.put(U) ; 

Text.Io ,new_line; 

Text_Io.put_line("Effective  height  of  the  plume  is  ") ; 
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My_Float_Io . put (H) ; 

Text_Io.new_line; 

Textile. put_line("Stability  condition  given  is  "); 
Text_Io. put (Stable) ; 

Text_Io.new_line; 

Text_Io .put_line("The  calculated  concentration  is  "); 
My _Float_Io .put (Concentration) ; 

Text_Io .new.line; 

end  Print. Info; 


—  PROCEDURE:  gausmod  (a.k.a  main) 

—  DESCRIPTION:  This  is  the  main  program  and  calls  all  the  above 

modules 

—  INPUT  PARAMETERS:  none 

—  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES:  none 

—  GLOBALS  USED:  All. 

—  CALLED  BY:  none 

—  CALLS:  Get.Info,  Calc.PG.Sigma,  Calc.Concentration, 

and  Print. Info. 


begin  — MAIN 

Get.Info(Q,  X,  Y,  Z,  U,  H,  Stable); 
Calc.PG.Sigma(Stable,  X,  SigmaY,  SigmaZ) ; 
Calc.Concentration(SigmaY ,  SigmaZ,  Q,  Y,  Z,  U,  H, 

Concentration) ; 

Print. Info (Concentration,  X,  Y,  Z,  Q,  U,  H,  Stable); 
end  gausmod; 
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Appendix  B.  Numerical  solutions  Ada  code 


The  B-series  appendices  contain  the  Ada  programming  language  programs 
which  contain  the  finite  difference  schemes  used  to  solve  the  advection  equation,  one-, 
two-,  and  three-dimensional  advection-diffusion  equations,  and  the  two-dimensional 
steady-  state  equation. 

Appendix  B.l  has  the  program  for  the  advection  equation.  Appendix  B.2 
has  the  program  for  the  1-D  advection-diffusion  equation.  Appendix  B.3  has  the 
program  for  the  2-D  advection-diffusion  equation.  Appendix  B.4  has  the  program 
for  the  3-D  advection-  diffusion  equation.  Appendix  B.5  has  the  program  for  the 
2-D  steady-state  equation. 

The  description  for  the  variables  used  in  the  programs  is  within  each  program 
in  the  form  of  comment  blocks  before  each  procedure  within  the  program.  Comment 
lines  start  with  a  double  hyphen  and  comment  blocks  begin  and  end  with  a  line  of 
hyphens. 
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B.l  Advection  Equation  Ada  Code 

—  FILE:  exl31.a 

—  PROJECT:  Thesis  work 

—  DATE:  15  Jul  93 

—  VERSION:  Version  1.0 

—  AUTHORS:  Capt  Dave  Paal 

—  DESCRIPTION:  This  program  solves  a  hyperbolic  partial 

differential  equation  using  the  Lzut-Friedrichs  scheme. 
The  partial  differential  equation  to  be  solved  is: 

d/dt  (U)  +  d/dx  (U)  =  0 

—  OPERATING  SYSTEM:  UNIX/ Sun  Sparc  Station 

—  LANGUAGE:  Meridian  Ada 

—  FILES  USED:  Output:  niiml31.dat. 


—  CONTEXT  CLAUSES 

with  text.io; 
with  Math.Lib; 
use  Math_Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  probl31  is 

—  TYPE  DECLARATION 


type  Vector  is  array  (1..50)  of  float 


—  GLOBAL  VARIABLES  AND  EXCEPTIONS 


V.Old 

Vector; 

V.New 

Vector; 

Delta_X 

float ; 

Delta_T 

float ; 

Min_X 

float ; 

Max_X 

float ; 

T_Value 

float; 

T 

float ; 

Outf ile 

Text_Io .File_Type; 

—  PROCEDURE:  Get  information  for  calculation 

—  DESCRIPTION:  This  procedure  gets  the  values  of  Delta.X, 

Delta_T,  Min_X,  Max_X,  and  Iterations. 

—  INPUT  PARAMETERS:  None. 

—  OUTPUT  PARAMETERS:  Del_X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 

Low_X  :  lower  bound  for  X 

Hi_X  :  upper  bound  for  X 

Time  :  time  value  wanted  for  calculation 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


procedure  Get_Info  (Del_X  :  in  out  float; 

Del_T  :  in  out  float; 

Low_X  :  in  out  float; 

Hi_X  :  in  out  float; 

Time  :  in  out  float)  is 

begin 

Text_Io.put_line("Enter  the  value  (Floating  point,  ie  0.1)  of 
Delta  X.’’); 

My_Float_Io . get (Del_X) ; 

Text_Io .put_line("Enter  the  value  (Floating  point,  ie  0.1)  of 
Delta  T."); 

My_Float_Io.get (Del_T) ; 

Text_Io.put_line( "Enter  the  minimum  value  (Floating  point,  ie 
-2.0)  of  X."); 

My_Float_Io.get(Low_X) ; 

Text.Io .put_line("Enter  the  maximum  value  (Floating  point,  ie 
3.0)  of  X."); 

My_Float_Io.get (Hi_X) ; 

Text_Io .put_line("Enter  calculation  time  value  (Floating  point, 
ie  1.6)."); 

My_Float_Io.get (Time) ; 
end  Get.Info; 


—  PROCEDURE:  Initial  Vector 

—  DESCRIPTION:  This  procedure  initializes  the  vector 
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using  Del_X,  Low_X,  Hi_X.  If  |x|<=1.0  the  vector  element 
is  given  the  value  l-|x|,  otherwise  the  vector  element  is 
zero . 

—  INPUT  PARAMETERS:  V_01d  :  vector  to  represent  x  direction 

Del_X  :  space  increment er 
Del_T  :  time  increment er 
Low_X  :  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 

—  OUTPUT  PARAMETERS:  V_01d  :  initial  vector  for  problem 

—  LOCAL  VARIABLES:  X  :  Value  of  X  for  each  iteration  of  loop 

Count  :  number  of  X  increments. 

—  GLOBALS  USED:  None. 

—  CALLED  BY :  main . 

—  CALLS :  none . 


procedure  Initial.Vector  (V_01d  :  in  out  Vector; 

Del_X,  Low_X,  Hi_X  :  in  float)  is 

X  :  float  :=  Low_X; 

count  :  integer  :=  integer((Hi_X  -  Low_X)/Del_X) ; 
begin 

for  i  in  1.. count  loop 
if  abs(X)  <=  1.0  then 
V.Old(i)  :=  1.0  -  abs(X); 
else 

V_01d(i)  :=  0.0; 
end  if; 

X  :=  X  +  Del.X; 
end  loop; 

end  Initial.Vector; 


—  PROCEDURE:  Compute  New  Vector 

—  DESCRIPTION:  This  procedure  calculates  the  solution  vector 

in  increments  each  time  it  is  called  by  the  main 
program . 

—  INPUT  PARAMETERS:  V_01d  :  Vector  for  calculating  amswer 

V_New  :  Vector  for  calculating  euxswer 
Del_X  :  incrementor  for  space 
Del_T  :  incrementor  for  time 
Low_X  :  lower  bound  for  X 
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Hi_X  :  upper  bound  for  X 

—  OUTPUT  PARAMETERS:  V_01d  :  Becomes  vector  with  values  for  this 

iteration 

V_Neu  :  Vector  for  calculating  answer 

—  LOCAL  VARIABLES:  Lzuabda  :  float  Delta_T/Delta_X 

Count  :  Counter  for  Vector  array 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 

procedure  Compute_New(V_01d  ;  in  out  Vector; 

V_New  :  in  out  Vector; 

Del_X,  Del_T,  Low_X,  Hi_X  :  in  float)  is 
Lambda  :  float  :=  Del_T/Del_X; 

Count  :  integer  :=  integer ((Hi_X  -  Low_X)/Del_X) ; 
begin 

V.New(l)  :=0.0; 

for  i  in  2.. (Count  -I)  loop 

V_New(i)  :=  0.5*(V.01d(I+l)  +  V.Old(I-i))  - 

0.5*Lambda*(V_01d(l+l)  -  V_01d(I-l)); 

end  loop; 

V_New(Count)  :=  V_New (Count-1 ) ; 
for  i  in  1 . . Count  loop 
V_01d(i)  :=  V.New(i); 
end  loop; 
end  Compute.New; 


—  PROCEDURE:  Print. Info 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 

—  It  prints  out  the  vector  with  the  points  for  the  time 

—  iteration  given. 

—  INPUT  PARAMETERS:  V.Old  :  Vector  for  calculating  answer 

V.New  :  Vector  for  calculating  einswer 
Del.X  :  incrementor  for  space 
Del.T  :  incrementor  for  time 
Low.X  :  lower  bound  for  X 
Hi.X  :  upper  bound  for  X 

—  OUTPUT  PARAMETERS:  V.New  :  Vector  euiswer 
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—  LOCAL  VARIABLES:  Count  :  integer  counter  for  array 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


procedure  Print. Info (V_New  ;  in  out  Vector; 

Del.X,  Del.T,  Low.X,  Hi.X,  Time  :  in  float; 

File  :  in  out  Text.Io.File.Type)  is 

Count  :  integer  :=  integer((Hi_X  -  Low.X) /Del.X) ; 

X.Value  :  float  :=  Low.X+Del.X; 
begin 

Text.Io.put.lineC'The  following  are  the  results  using  the 
Lax-  "); 

Text. lo. put. lineC'Friedrichs  scheme  on  a  hyperbolic  partial 
difeq.") ; 

Text.Io.put.lineC'The  resulting  vector  is  "); 
Text.Io.new.line(File) ; 

Text.Io.put_line(File,"#  The  following  are  the  results  using 
the  ") ; 

Text.Io.put_line(File,"#  Lax-Friedrichs  scheme  on  a 
hyperbolic") ; 

Text.Io.put_line(File,"#  partial  difeq.  The  resulting  vector 
is  ") ; 

Text. lo. put (File,"#  Output  for  Min.x*") ; 

My. Float. lo. put (File, Low.X, 3, 2,0) ; 

Text.Io.new.line(File) ; 

Text. lo  .put  (File,  "#  Mauc.x=")  ; 

My.Float.Io.put (File, Hi.X, 3, 2,0) ; 

Text.Io .new.line(File) ; 

Text. lo. put (File,"#  Delta.X=") ; 

My.Float.Io.put (File, Del.X, 3, 5,0) ; 

Text.Io.new.line(File) ; 

Text.Io .put(File,"#  Delta.T=") ; 

My . Float. lo .put (File, Del.T, 3, 5,0) ; 

Text.Io .new.line(File) ; 

Text.Io .put (File,"#  Time  of  calculation=") ; 

My. Float. lo . put (File , Time ,3,5,0); 

Text.Io .new.line(File) ; 
for  i  in  1.. Count  loop 

My.Float.Io.put(X.Value, 1,2,0) ; 

Text.Io .put ("  ")  ; 
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'Iy.Float_Io.put(V_New(i)  ,1,6,2) ; 
Text_Io.new_line; 
Text_Io.8et_col(File,l) ; 

My_Float_Io .put (File,X_Value, 1 ,2,0); 
Text_Io.set_col(File,10) ; 

My _Float_Io.put(File,V_Mew(i) ,1,6,2) ; 
X_Value  :=  X.Value  +  Del.X; 
end  loop; 

Text.Io.new.line; 

Text_Io .new_line; 
end  Print.Info; 


—  PROCEDURE:  Main 

—  DESCRIPTION:  This  is  the  main  program  and  calls  all  the  above 

modules 

—  INPUT  PARAMETERS:  none 

—  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES;  none 

—  GLOBALS  USED:  All. 

—  CALLED  BY:  none 

—  CALLS:  Get.Info,  Initial.Vector,  Compute.New,  Print.Info. 


begin  —MAIN 

Text.Io . Create (Outf ile .Text _Io . Out .File , "numl31 . dat " ) ; 
Get.Info  (Delta.X,  Delta.T,  Min_X,  Meuc.X,  T.Value) ; 
Initial.Vector  (V.Old,  Delta.X,  Min.X,  Meuc.X) ; 

T  :=  Delta.T; 

while  T  <=  T.Value  loop 

Compute.NewC V.Old,  V.New,  Delta.X,  Delta.T,  Min.X,  Max.X) ; 
T  :=  T  +  Delta.T; 
end  loop; 

Print.Info (V.New,  Delta.X,  Delta.T,  Min.X,  Meuc.X,  T.Value, 
Outf ile)  ; 

Text. Io.Close(Outf ile) ; 
end  problSl; 
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B.2  1-D  Advection-Diffusion  Equation  Ada  Code 

—  FILE:  exadvdif.a 

—  PROJECT:  Thesis  work 

—  DATE:  19  Jul  93 

—  VERSION:  Version  1.0 

—  AUTHORS:  Capt  Dave  Paal 

—  DESCRIPTION:  This  program  solves  a  hyperbolic  partial 

differential  equation  using  the  forward-time  auid 
central-space  forward  difference  scheme. 

The  partial  differential  equation  to  be  solved  it: 

d/dt  (U)  +  d/dx  (U)  =  b  ♦  d**2/dx**2  (U) 

"  OPERATING  SYSTEM:  UNIX/Sun  Sparc  Station 

—  LANGUAGE:  Meridian  Ada 

—  FILES  USED:  Output:  advdifl.dat. 


—  CONTEXT  CLAUSES 

with  text.io; 
with  Math.Lib; 
use  Math.Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  exadvl  is 

—  TYPE  DECLARATION 

type  Vector  is  array  (1..50)  of  float; 


—  GLOBAL  VARIABLES  AND  EXCEPTIONS 


V.Old 

Vector; 

V_New 

Vector; 

Delta_X 

float ; 

Delta_T 

float ; 

Min_X 

float ; 

Max.X 

float ; 

T_Value 

float ; 

A,  B 

float ; 
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Bmu  :  float; 

T  :  float ; 

Outfile  :  Text_Io .File.Type; 

—  PROCEDURE:  Get  information  for  calculation 

—  DESCRIPTION:  This  procedure  gets  the  values  of  Del_X,  Del_T, 

Low_X,  Hi_X,  jmd  Time. 

—  INPUT  PARAMETERS:  None. 

—  OUTPUT  PARAMETERS:  Del_T  :  time  incrementor 

Del_X  :  space  incrementor 
Low_X  :  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 
A,  B  :  constants  (positive) 

Time  :  time  calculation  is  wauited 

—  LOCAL  VARIABLES:  None. 

"  GLOBALS  USED: 

—  CALLED  BY:  main. 

—  CALLS :  none . 

procedure  Get. Info  (Del.X  :  in  out  float; 

Del_T  :  in  out  float; 

Low.X  :  in  out  float; 

Hi.X  :  in  out  float; 

A,  B  :  in  out  float; 

Time  :  in  out  float)  is 

begin 

Text.Io .put_line("Enter  the  value  (Floating  point,  ie  0.1)  of 
Delta  X."); 

My_Float_Io.get(Del_X) ; 

Text  _Io.  put  .line  ("Enter  the  value  (Floatingpoint,  ie  0.05)  of 
Delta  T."); 

My.Float_Io.get (Del.T) ; 

Text _Io.put_line( "Enter  the  minimum  value  (Floating  point,  ie 
-2.0)  of  X."); 

My_Float.Io.get(Low_X) ; 

Text.Io.put.lineC'Enter  the  meiximum  value  (Floating  point,  ie 
3.0)  of  X."); 

My.Float.Io.get (Hi.X) ; 

Text.Io.put.lineC'Enter  the  value  (Positive  Floating  point,  ie 
0.5)  of  A."); 


My_Float_Io.g0t (A) ; 

Text_Io.put_line("Enter  the  value  (Positive  Floating  point,  ie 
0.5)  of  B."); 

My_Float_Io.get (B) ; 

Text_Io.put_line("Enter  calculation  time  value  (Floating  point, 
ie  1.6)."); 

My.Float.Io. get (Time) ; 
end  Get.Info; 


—  PROCEDURE:  Initial  Vector 

—  DESCRIPTION:  This  procedure  initializes  the  vector 

using  Del_X,  Low.X,  Hi_X.  If  lxl<=1.0  the  vector  element 
is  given  the  value  l-|x|,  otherwise  the  vector  element  is 
zero . 

—  INPUT  PARAMETERS:  Del_X  :  space  incrementor 

Low_X,  Hi_X  :  min  and  max  range 
A,  B  :  positive  consteoits 

—  OUTPUT  PARAMETERS:  V_01d  :  V(0) 

—  LOCAL  VARIABLES:  X  :  float 

Count  :  number  of  X  increments 
Mu  :  Del_T/(Del_X  ♦  Del.X) 

—  GLOBALS  USED:  none. 

—  CALLED  BY:  main. 

--  CALLS :  none . 


procedure  Initial.Vector  (V_01d  :  in  out  Vector; 

Del_X,  Del.T,  Low_X,  Hi_X,  A,  B  :  in  float)  is 

X  :  float  :=  Low_X; 

count  :  integer  :=  integer((Hi_X  -  Low_X)/Del_X) ; 

Mu  :  float  :=  Del_T/(Del_X  ♦  Del_X) ; 

begin 

if  B*Mu  >  0.5  then 

Text_Io.put_line("B*Mu  >  0.5  causes  unstable  results"); 
else 

for  i  in  l..co\int  loop 
if  abs(X)  <=  1.0  then 

V_01d(i)  :=  1.0  -  abs(X); 
else 

V_01d(i)  :=  0.0; 
end  if; 


X  :=  X  +  Del_X; 
end  loop; 
end  if; 

end  Initial.Vector; 


—  PROCEDURE:  Compute  Vector  answer  for  time  wanted. 

—  DESCRIPTION:  This  procedure  calculates  new  vector 

—  INPUT  PARAMETERS:  V.Old  :  Vector  for  V(N) 


V_New  :  Vector  for  V(N+1) 
Del_X  :  incrementor  for  space 
Del_T  :  incrementor  for  time 
Low_X  :  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 
A,  B  :  positive  constants 
OUTPUT  PARAMETERS:  V.Old  :  Vector  for  V(N) 

V.New  :  Vector  for  V(N+1) 
LOCAL  VARIABLES:  Mu  :  float  Del_T/Del.X**2 

Count  :  Counter  for  Vector  array 
Alpha  :  (Del.X  ♦  A)/(2*B) 

GLOBALS  USED:  None. 

CALLED  BY:  main. 

CALLS :  none . 


—  LOCAL  VARIABLES: 


mam. 


procedure  Compute.New (V.Old  :  in  out  Vector; 

V.New  :  in  out  Vector; 

Del.X,  Del.T,  Low.X,  Hi.X,  A,  B  :  in  float)  is 

Count  :  integer  :=  integer((Hi.X  -  Low.X) /Del.X) ; 

Mu  :  float  :=  Del.T/ (Del.X  *  Del.X); 

Alpha  :  float  :=  (Del.X  *  A)/(2.0*B); 

begin 

V.New(l)  :=0.0; 

for  i  in  2.. (Count  -  1)  loop 

V.New(i)  :=  ((1.0  -  2.0  *  B  *  Mu)  ♦  V_01d(i))  +  (B  *  Mu  ♦ 

(1.0  -  Alpha)  *  V_01d(i+1))  + 

(B  *  Mu  *  (1.0  +  Alpha)  *  V.Old(i-l)); 

end  loop; 

V.New(Count)  :=  V.New (Count -1 ) ; 
for  i  in  1.. Count  loop 


V_01d(i)  :=  V_New(i); 
end  loop: 
end  Cofflpute_New; 


”  PROCEDURE:  Print.Info 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 

It  prints  the  point  and  corresponding  concentration  for  the 

—  wanted  time  to  both  the  screen  and  an  output  file  named 

—  advdifl.dat. 

—  INPUT  PARAMETERS:  V_New  :  Vector  answer 

Del.X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 

Low_X  :  lower  bound  for  X 

Hi_X  :  upper  bound  for  X 

—  OUTPUT  PARAMETERS:  V.New  :  Vector  answer 

—  LOCAL  VARIABLES:  Count  :  integer  counter  for  array 

X_Value  :  value  of  x  for  each  iteration 

through  vector 

Alpha  :  (Del.X  ♦A)  /  (2.0  ♦  B)  <=  1.0  for 

—  no  oscillation 

Bmu  :  (B  *  Del .T)/(Pal_X  +  Del.X)  <=0.5 

for  stability 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


procedure  Print.Info (V.New  :  in  out  Vector; 

Del.X,  Del.T,  Low.X,  Hi.X,  Time,  A,  B  :  in  float; 
File  :  in  out  Text.Io.File.Type)  is 

Count  :  integer  :=  integer ((Hi.X  -  Low.X) /Del.X) ; 

X.Value  :  float  :=  Low.X+Del.X; 

Alpha  :  float  :=  (Del.X  ♦A)  /  (2.0  *  B); 

Bmu  :  float  :=  (B  *  Del.T) /(Del.X  *  Del.X); 

begin 

Text. I o  put .line ("The  following  are  the  results  using  forward 
time  ") ; 

Text. I o. put .line ("central  space  scheme  on  a  hyperbolic  partial 

difeq. ") ; 

Text. lo .put. lineC'The  resulting  vector  is  ") ; 
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Text_Io.put_line(File,"#  The  following  are  the  results  using 


the  ") ; 

Text_Io.put_line(File,"#  forward  time/central  space  scheme  on 

a  *') ; 

Text_Io.put_line(File, "#  Ct  +  aCx  =  bCxx."); 
Text_Io.put_line(File,"#  The  resulting  vector  is:  ") ; 

Text_Io. put (File,"#  Output  for  Min_x="); 

My_Float_Io .put(File,Low_X,3,2,0) ; 

Text_Io.new_line(File) ; 

Text_Io .put (File, "#  Max_x=")  ; 

My_Float_Io .put(File,Hi_X,3,2,0) ; 

Text_Io.new_line(File) ; 

Text_Io. put (File,"#  Delta  X=") ; 

My_Float.Io .put(File,Del_A,3,5 ,0) ; 

Text_Io.new_line(File) ; 

Text.Io . put (File , "#  Delta_T=") ; 

My.Float^Io .put(File,Del_T,3,5,0) ; 

Text.Io .new_line(File) ; 

Text_Io. put (File,"#  Time  of  calculation=") ; 

My_Float_Io .put(File,Time,3,5  f) ; 

Text_Io.new_line(File) ; 

Text.Io. put (File,"#  Alpha®  (Dejta.X  ♦  A)/2B®"); 

My.Float_Io .put (File, Alpha, 3, 5,0) ; 

Text_Io.new_line(File) ; 

Text_Io. put (File,"#  Bmu=  B*Delta_T/Delta_X**2=") ; 

My_Float_Io .put(File,Bmu,3,5,0) ; 

Text.Io .new_line(File) ; 
for  i  in  1.. Count  loop 

Text_Io.set_col(File,l) ; 

My_Float_Io .put(File,X_Value,2,2,0) ; 

Text_Io . set_col(File, 15) ; 

My_Float_Io .put(File,V_New(i) ,2,6,2) ; 

X.Value  :=  X.Value  +  Del_X; 
end  loop; 

T  ext  _ I o . new_l ine ; 

T  ext  _ I o . new_l ine ; 
end  Print.Info; 


—  PROCEDURE:  Examplel.3.1  (a.k.a  main) 

—  DESCRIPTION:  This  is  the  main  prograun  auid  calls  all  the  above 
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modules 

—  INPUT  PARAMETERS:  none 

—  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES:  none 

—  GLOBALS  USED:  All. 

—  CALLED  BY:  none 

—  CALLS:  Get_Info,  Initial .Vector,  Compute.New,  Print.Info. 


begin  --MAIN 

Text.Io . Create (Outfile ,Text_Io . Out.File , "advdif 1 . dat ") ; 
Get.Info  (Delta.X,  Delta.T,  Min_X,  Mauc.X,  A,  B,  T.Value) ; 
Initial.Vector  (V_01d,  Delta.X,  Delta.T,  Min_X,  Max_X,  A,  B) ; 
Bmu  :=  B*Delta_T/(Delta_X*Delta_X) ; 
if  Bmu  >  0.5  then 

Text.Io .put_line("Unstable  result,  try  again."); 
else 

T  :=  Delta.T; 

while  T  <=  T.Value  loop 

Compute_New(V.Dld,  V.New,  Delta.X,  Delta.T,  Min_X,  Max.X, 

A,  B); 

T  :=  T  +  Delta.T; 
end  loop; 

Print_Info(V_New,  Delta.X,  Delta.T,  Min_X,  Max.X,  T.Value, 
A,  B,  Outfile); 

end  if; 

Text_Io.Close(Outfile) ; 
end  exadvl; 
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B.S  2-D  Advection- Diffusion  Equation  Ada  Code 

—  FILE:  exadvdif2.a 

—  PROJECT:  Thesis  work,  example  calculation 

—  DATE:  23  Jul  93 

—  VERSION:  Version  1.0 

—  AUTHORS:  Capt  Dave  Paal 

—  DESCRIPTION:  This  progreim  solves  a  partial  differential 

equation  using  the  forward  time  and  central  space. 
The  partial  differential  equation  to  be  solved  is: 

Ct  +  A*Cx  =  Kl*Cxx  +  K2*Cyy 

—  OPERATING  SYSTEM:  UNIX/Sun  Sparc  Station 

—  LANGUAGE:  Meridian  Ada 

—  FILES  USED:  Output:  expltest.dat 


—  CONTEXT  CLAUSES 

with  text.io; 
with  Math_Lib; 
use  Math_Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  ctprob  is 

—  TYPE  DECLARATION 

type  Matrix  is  array  (1 . .80,1 . .80)  of  float; 


—  GLOBAL  VARIABLES  AND  EXCEPTIONS 


C_N 

Matrix; 

C.New 

Matrix; 

Delta_X 

float ; 

Delta.Y 

float ; 

Delta_T 

float ; 

Min.X 

float ; 

Max.X 

float ; 

Min_Y 

float ; 

Max_Y 

float ; 

A 

:  float; 

Kl 

:  float; 

K2 

:  float; 

T.Value 

:  float; 

T 

:  float; 

Outf ile 

;  Text_Io .File_Type; 

—  PROCEDURE:  Get  information  for  calculation 

—  DESCRIPTION:  This  procedure  gets  the  values  of  Del_X,  Del_T, 

Del.Y,  Lotf.X,  Hi_X,  Low.Y,  Hi_Y,  A,  Kl,  K2,  and  Time. 

”  INPUT  PARAMETERS:  None. 

—  OUTPUT  PARAMETERS:  Del_X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 
Del_Y  :  incrementor  for  Y  direction 
Lotf_X  :  lover  bound  for  X 
Hi_X  :  upper  bound  for  X 
Lov_Y  :  lower  bound  for  Y 
Hi_Y  :  upper  bound  for  Y 
A  :  constant  wind  sj  -ed 

Kl,  K2  :  diffusivity  constants 
Time  :  wanted  time  level 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED:  Del.X,  Del.T,  Del.Y,  Low.X,  Hi_X,  Low.Y,  Hi.Y, 

A,  Kl,  K2,  and  Time. 

—  CALLED  BY :  main . 

—  CALLS :  none . 

procedure  Get_Info  (Del_X  :  in  out  float; 

Del_Y  :  in  out  float; 

Del_T  :  in  out  float; 

Low_X  :  in  out  float; 

Hi_X  :  in  out  float; 

Low_Y  :  in  out  float; 

Hi_Y  :  in  out  float; 

A  :  in  out  float; 

Kl,  K2:  in  out  float; 

Time  :  in  out  float)  is 

begin 

Text_Io .put_line("Enter  the  value  (Floating  point)  of  Delta 
X.  "); 


My_Float_Io.get(Del_X) ; 

Text_Io.put_line("Eiiter  the  value  (Floating  point)  of  Delta 
Y.  "); 

My_Float_Io.get(Del_Y) ; 

Text_Io.put_line("Enter  the  value  (Floating  point)  of  Delta 
T.  "); 

My_Float_Io.get(Del_T) ; 

Text_Io.put_line("Enter  the  minimum  value  (Floating  point) 
of  X.  "); 

My_Float_Io.get(Low_X) ; 

Text_Io.put_line("Enter  the  maximum  value  (Floating  point) 
of  X.  "); 

My  .Float _Io. get (Hi_X) ; 

Text _ I 0. put .line ("Enter  the  minimum  value  (Floating  point) 
of  Y.  "); 

My.Float_Io.get(Low.Y) ; 

Text. lo. put. line("Enter  the  maximum  value  (Floating  point) 
of  Y.  "); 

My.Float.Io .get (Hi.Y) ; 

Text. lo. put. line("Enter  the  value  (Positive  Floating  point) 
of  A.  ") ; 

My.Float.Io. get (A) ; 

Text. lo. put. line("Enter  the  value  (Positive  Floating  point) 
of  Kl.  "); 

My.Float.Io .get (Kl) ; 

Text. lo .put. lineC'Enter  the  value  (Positive  Floating  point) 
of  K2.  "); 

My.Float.Io .get (K2) ; 

Text. lo. put. lineC'Enter  the  number  of  time  value  wanted 
(float) .") ; 

My.Float.Io. get (Time) ; 
end  Get. Info; 


—  PROCEDURE:  Initial  Matrix 

—  DESCRIPTION:  This  procedure  initializes  the  matrix 

using  Del.X,  Del.Y,  Low.X,  Hi.X,  Low.Y,  Hi.Y. 

—  INPUT  PARAMETERS:  C.N,  Del.X,  Del.Y,  Low.X,  Hi.X,  Low.Y,  Hi.Y 

—  OUTPUT  PARAMETERS:  C.N. 

—  LOCAL  VARIABLES:  X  :  X  value  incrementer. 

Y  :  Y  value  incrementer. 

Count. X  :  number  of  X  increments. 

Count. Y  :  number  of  Y  increments. 
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—  PROCEDURE:  Compute  final  Matrix 

—  DESCRIPTION:  This  procedure  calculates  new  vector 

—  INPUT  PARAMETERS:  C_N  :  Matrix  for  C_N(i,j) 

C_New  :  Matrix  for  calculating  answer 

Del_X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 

Del_Y  :  incrementor  for  Y  direction 

Low_X  :  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 
Low_Y  :  lower  bound  for  Y 
Hi_Y  :  upper  bound  for  Y 


A  :  constant 

K1  :  X  diffusivity  constant 

K2  :  Y  diffusivity  constant 

—  OUTPUT  PARAMETERS:  C_N  :  Matrix  for  C(N) 

C_New  :  Matrix  for  calculating  answer 

—  LOCAL  VARIABLES:  Mu  :  (A*DQl_T)/(2.0*Del_X) 

Alpha  :  (Kl*Del_T)/(Del_X*Del_X) 

Beta  :  (K2*Del.T)/(Del.Y*Del_Y) 

Count_Y  :  integer ( (Hi_Y  -Low_Y)/Del_Y) 
Count _X  :  integer ((Hi_X  -  Low_X)/Del_X) 

—  GLOBALS  USED:  C.N,  C_New.  Del.X,  Del_T,  Del_Y,  Lou.X,  Hi_X 

Low_Y.  Hi_Y,  A,  Kl,  K2. 

—  CALLED  BY :  main . 

—  CALLS :  none . 

"  ANALYSIS:  0(1). 


procedure  Compute_Final(C_N  :  in  out  Matrix; 

C.New  :  in  out  Matrix; 

Del.X,  Del.T,  Del.Y  :  in  float; 

Low.X,  Hi.X,  Low.Y,  Hi.Y  :  in  float; 

A,  Kl,  K2  :  in  float)  is 

Mu  :  float  :*  (A*Del.T)/(2.0*Del.X) ; 

Alpha  :  float  :=  (Kl*Del.T)/(Del.X*Del_X) ; 

Beta  :  float  :=  (K2*Del_T)/(Del.Y*Del_Y) ; 

Count _X  :  integer  :=  integer ((Hi.X  -  Low.X) /Del.X) ; 

Count _Y  :  integer  :=  integer((Hi_Y  -  Low.Y) /Del.Y) ; 

begin 

for  i  in  l..Count.X  loop 
C.New(i ,Co\int_Y)  :=0.0; 

C.New(i,l)  :=0.0; 
end  loop; 

for  i  in  l..Count_Y  loop 
C.New (Count. X,i)  :=0.0; 

C.New(l,i)  :=0.0; 
end  loop; 

for  j  in  2..(Count.Y  -  1)  loop 
for  i  in  2.. (Count .X  -  1)  loop 

C.New(i,j)  :=  C.N(i,j)  -  Mu*(C_N(i+l, j)  -  C_N(i-l,j))  + 
Alpha* (C.N (i+l,j)  -  2.0*C_N(i,j)  + 
C_N(i-l,j))  + 
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Beta*(C_N(i , j+1)  -  2.0*C_N(i,j)  + 

c_N(i,j-i)); 

end  loop; 
end  loop; 

for  i  in  l,.Count_X  loop 
for  j  in  l..Count_Y  loop 
C_N(i, j)  :=  C_New(i, j)  ; 
end  loop; 
end  loop; 

end  Compute_Final ; 


”  PROCEDURE:  Print.Info 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 

—  It  prints  out  the  matrix  with  the  points  for  the  time 

—  iteration  given. 

—  INPUT  PARAMETERS:  C_New  :  Matrix  for  calculating  auiswer 

Del_X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 

Del_Y  :  incrementor  for  Y  direction 

Low_X  :  lower  bound  for  X 

Hi_X  :  upper  bound  for  X 

Low_Y  :  lower  bound  for  Y 

Hi_Y  :  upper  bound  for  Y 

—  OUTPUT  PARAMETERS:  C.New  :  Matrix  answer 

File  :  output  file 

—  LOCAL  VARIABLES:  Count_X,  Count _Y  :  matrix  integer  counters 

X.Value,  Y_Value  :  X  and  Y  points 

—  GLOBALS  USED:  C.New,  Del_X,  Del.T,  Del_Y,  Low.X,  Hi.X, 

Low_Y,  Hi_Y,  File. 

—  CALLED  BY :  main . 

—  CALLS :  none . 


procedure  Print.Info (C.New  :  in  out  Matrix; 

Del.X,  Del.T,  Del.Y  :  in  float; 

Low.X,  Hi.X,  Low.Y,  Hi.Y  :  in  float; 
File  :  in  out  Text.Io.File.Type)  is 

Count.X  :  integer  :=  integer((Hi.X  -  Low.X) /Del.X) ; 
Count.Y  :  integer  :=  integer((Hi_Y  -  Low.Y) /Del.Y) ; 

X. Value  :  float  :=  Low.X; 

Y. Value  :  float  :=  Low.Y; 
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begin 

Text_Io.put_line("iThe  following  are  the  results  using  the 
forward  time  "); 

Text_Io.put_line("#  central  space  for  the  partial 
diferential  equation:"); 

Text_Io.put_line("#  Ct  +  A*Cx  =  Kl*Cxx  +  K2*Cyy  ") ; 

Text_Io.put_line("#  The  resulting  matrix  is:  ") ; 

Text.Io . new_line ; 

Text.Io .new_line(File,2) ; 

Text_Io.put_line(File,"#  Output  for  a=.5,  kl=.5,  k2=.5"); 
Text_Io.put_line(File,"#  delta_x=.2,  delta_y=.2, 
delta_t=.02,t=.06") ; 

for  j  in  1.. Count _Y  loop 
X_Value  :=  Low_X; 
for  i  in  l..Count_X  loop 
Text_Io.set_col(File,l) ; 

My_Float_Io.put (File.X.Vadue, 1 ,2,0) ; 
Text_Io.set_col(File,10) ; 

My.Float_Io. put (File,  Y.Value, 1,2,0) ; 
Text_Io.set_col(File,20) ; 
My_Float_Io.put(File,C_Mew(i. j) .1.6»2) ; 
Text_Io.new_line(File) ; 

X.Value  :=  X_Value  +  Del.X; 
end  loop; 

Y. Value  :=  Y_ Value  +  Del.Y; 
end  loop; 

Text_Io.new_line; 

Text.Io .new.line; 
end  Print.Info; 


—  PROCEDURE:  Main 

—  DESCRIPTION:  This  is  the  main  prograun  and  calls  all  the  above 

modules 

—  INPUT  PARAMETERS:  none 

—  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES:  none 

—  GLOBALS  USED:  All. 

—  CALLED  BY:  none 

—  CALLS:  Get.Info,  Initial_Matrix,  Compute_Final ,  Print.Info. 


begin  — MAIN 
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Text  _ I o . Great  e(0utfile, Text  _ I o . Out  _F ile , " explt  est . dat " ) ; 
Get_Info  (Delta.X,  Delta.Y,  Delta.T,  Min.X,  Max.X,  Min_Y, 
Max.Y.  A,  Kl,  K2,  T_ Value) ; 

Initial_Matrix  (C_N,  Delta_X,  Delta.Y,  Min_X,  Max^X,  Min_Y , 
Max.Y) ; 

T  :=  Delta.T; 

while  T  <=  T.Value  loop 

Coinpute.Final(C.M,  C.New,  Delta.X,  Delta.T,  Delta.Y,  Min.X, 
Max.X,  Min.Y,  Max.Y,  A,  Kl,  K2) ; 

T  :=  T  +  Delta.T; 
end  loop; 

Print.Info (C.New,  Delta.X,  Delta.T,  Delta.Y,  Min.X,  Max.X, 
Min.Y,  Max.Y,  Outfile); 

Text. lo.CloseCOutf ile) ; 
end  ctprob; 
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B.4  3-D  Advection- Diffusion  Equation  Ada  Code 

—  FILE:  threedprob . a 

—  PROJECT:  Thesis  work 

—  DATE:  12  Aug  93 

—  VERSION:  Version  1.0 

—  AUTHORS:  Capt  Dave  Paal 

—  DESCRIPTION:  This  progreun  solves  a  partial  differential 

equation  using  the  forward  time  and  central  space. 

The  partial  differential  equation  to  be  solved  is: 

Ct  +  A*Cx  =  Kl*Cxx  +  K2*Cyy  +  K3*C2z  +  F(X,Y,Z) 

Where  F(X,Y,Z)  is  the  source  term. 

For  this  example  the  following  conditions  must  hold  for  a 
stable  non-oscillating  solution: 

(Ki*Del_T)/(Del_X*Del_X)<=  1/8 
and  (A*Del_X)/(2*Ki)<=  1.0 

where  the  i  subscript  is  1,2, or  3  for  X,Y,or  Z 
diffusivity  terms. 

These  conditions  limit  the  usefulness  of  this  program 
because  they  more  or  less  make  the  diffusivity  terms 
constant  when  in  reality  they  are  variable 
(K(i)»sigma(i)**2/(2*time) . 

Another  limitation  of  this  program  is  that  the  source  term 
is  hard  coded  as  F(0,0,Stack_Height)  and  the  program  is 
dependent  on  whether  the  Delta^X,(Y  and  Z)  terms  and  the 
minimum  values  of  each  will  increment  to  X=0,  Y=0, 
euid  Z=  height. 

—  OPERATING  SYSTEM:  UNIX/Sun  Sparc  Station 

—  LANGUAGE:  Meridian  Ada 

—  FILES  USED:  Output:  test3d.dat 


—  CONTEXT  CLAUSES 

with  text.io; 
with  Math.Lib; 
use  Math_Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  threed  is 
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TYPE  DECLARA'^ION 


type  Tri_Matrix  is  array  (0. .50,0. .50,0. .50)  of  float; 


—  GLOBAL  VARIABLES  AND  EXCEPTIONS 


C_N 

Tri.Matrix; 

C.New 

Tri.Matrix; 

C_Calc 

Tri.Matrix; 

Delta_X 

float 

Delta. Y 

float 

Delta.Z 

float 

Delta.T 

float 

Min.X 

float 

Meix.X 

float 

Min.Y 

float 

Max_Y 

float, 

Min.Z 

float 

Max.Z 

float 

A 

float 

T.Value 

float 

Q.Value 

float 

Stk.Hgt 

float 

T 

float 

Outf ile 

Text.Io .File.Type; 

—  PROCEDURE:  Get  information  for  calculation 

—  DESCRIPTION:  This  procedure  gets  the  values  of  Del_X,  Del_T, 

Del_Y,  Del_Z,  Low_X,  Hi_X,  Low_Y,  Hi_Y,  Low_Z,  Hi_Z, 

A,  euid  Time. 


—  INPUT  PARAMETERS:  None. 

—  OUTPUT  PARAMETERS:  Del_X 

Del_T 

Del_Y 

Del.Z 

Low_X 

Hi_X 

Low.Y 

Hi.Y 

Low_Z 


:  incrementor  for  space 
:  incrementor  for  time 
:  incrementor  for  Y  direction 
:  incrementor  for  Z  direction 
:  lower  bound  for  X 
:  upper  bound  for  X 
:  lower  bound  for  Y 
:  upper  bound  for  Y 
:  lower  bound  for  Z 
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Hi_Z  :  upper  bound  for  Z 
A  ;  const2mt  wind  speed 

Time  :  wanted  time  level 

—  LOCAL  VARIABLES:  None. 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


procedure  Get.Info  (Del_X  :  in  out  float; 

Del_Y  :  in  out  float; 

Del_Z  :  in  out  float; 

Del_T  :  in  out  float; 

Low_X  :  in  out  float; 

Hi_X  :  in  out  float; 

Low_Y  :  in  out  float; 

Hi_Y  :  in  out  float; 

Low_Z  :  in  out  float; 

Hi_Z  :  in  out  float; 

A  :  in  out  float; 

Time  :  in  out  float; 

Quantity  :  in  out  float; 

Height  :  in  out  float)  is 

begin 

Text_Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

X.  "); 

My_Float_Io . get (Del_X) ; 

Text.Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

Y.  "); 

My„Float_Io.get(Del_Y) ; 

Text_Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

Z.  "); 

My_Float_Io .get (Del_Z) ; 

Text.Io .put_line("Enter  the  value  (Floating  point)  of  Delta 
T.  "); 

My_Float_Io . get (Del_T) ; 

Text_Io.put_line("Enter  the  minimum  value  (Floating  point) 
of  X.  "); 

My_Float_Io .get (Low_X) ; 

Text.Io .put_line("Enter  the  maximum  value  (Floating  point) 
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of  X.  "); 

My_Float_Io .get(Hi_X) ; 

Text_Io .put_line("Enter  the  minimum  value  (Floating  point) 
of  Y.  "); 

My_Float_Io .get (Low_Y) ; 

Text_Io .put_line("Enter  the  majcimum  value  (Floating  point) 
of  Y.  "); 

My_Float_Io . get (Hi_Y) ; 

Text_Io.put_line( "Enter  the  minimum  value  (Floating  point) 
of  Z.  "); 

My_Float_Io .get (Low_Z) ; 

Text.Io .put_line("Enter  the  maximum  value  (Floating  point) 
of  Z.  "); 

My_Float_Io.get(Hi_Z) ; 

Text_Io .put_line("Enter  the  value  (Positive  Floating  point) 
of  A.  "); 

My_Float_Io . get (A) ; 

Text.Io .put_line("Enter  the  number  of  time  value  wanted 
(float)."); 

My_Float_Io .get (Time) ; 

Text_Io.put_line("Enter  the  value  (Positive  Floating  point) 
of  Q.  "); 

My _Float_Io. get (Quantity) ; 

Text_Io.put_line("Enter  the  height  of  the  stack  (float)."); 

My_Float_Io. get (Height) ; 
end  Get _ Info; 


—  PROCEDURE:  Compute  final  Tri.Matrix 

—  DESCRIPTION:  This  procedure  calculates  final  answer  through 

iterations  using  a  finite  difference  scheme  with  forward 
time/central  space. 

—  INPUT  PARAMETERS:  C_N  :  Tri.Matrix  for  C_N(i,j) 

C_New  :  Tri_Matrix  for  calculating  answer 

Del_X  :  incrementor  for  space 

Del_T  :  incrementor  for  time 

Del_Y  :  incrementor  for  Y  direction 

Del_Z  :  incrementor  for  Z  direction 

Low_X  :  lower  bound  for  X 

Hi_X  :  upper  bound  for  X 

Low_Y  :  lower  bound  for  Y 

Hi_Y  :  upper  bound  for  Y 
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Loh_Z 

Hi_Z 

A 

Time 

—  OUTPUT  PARAMETERS:  C_N 

C_New 

—  LOCAL  VARIABLES:  Mu 

K1 

K2 

K3 

Alpha 

Beta 

Gamma 

Count_Y 

Count _X 

Count _Z 

X 

Y 

Z 

Fofxyz 

Pi 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


lower  bound  for  Z 
upper  bound  for  Z 
consteint 

counter  for  time  iteration 
Tri_Matrix  for  C(N) 

Tri.Matrix  for  calculating  answer 
(A*Del_T) / (2 . 0*Del_X) 

X  diffusivity  constemt 
Y  diffusivity  constant 
Z  diffusivity  constant 
(Kl*Del_T)/(Del_X*Del_X) 
(K2*Del_T)/(Del_Y*Del_Y) 
(K3*Del_T)/(Del_Z*Del_Z) 
integer((Hi_Y  -Low_Y)/Del_Y) 
integer ((Hi_X  -  Low_X)/Del_X) 
integer((Hi_Z  -  Low_Z)/Del_Z) 
counter  for  X  axis 
counter  for  Y  axis 
counter  for  Z  2ucis 
source  forcing  function 
the  number  Pi 


procedure  Compute_Final(C_N  :  in  out  Tri.Matrix; 

C_New  :  in  out  Tri_Matrix; 

C_Calc  :  in  out  Tri_Matrix; 

T,  Del.X,  Del_Y,  Del.Z,  Del_T  :  in  float; 
Low_X,  Hi_X,  Low_Y  :  in  float; 

Hi_Y,  Low_Z,  Hi_Z  :  in  float; 

A,  Quauitity,  Height  :  in  float)  is 


K1  :  float  :=  1.0; 

K2  :  float  :=  1.0; 

K3  :  float  :=  1.0; 

Mu  :  float  :=  (A*Del_T)/(2.0*Del_X) ; 

Alpha  :  float  :=  (Kl*Del_T)/(Del_X*Del_X) ; 

Beta  :  float  :=  (K2*Del_T)/(Del_Y*Del_Y) ; 

Gamma  :  float  :=  (K3*Del_T)/(Del_Z*Del_Z) ; 

Count.X  :  integer  :=  integer((Hi_X  -  Low_X)/Del_X) ; 
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Count _Y  :  integer  :=  integer ((Hi_Y  -  Low_Y)/Del_Y) ; 
Count_Z  :  integer  :=  integer((Hi_Z  -  Low_Z)/Del_Z) ; 
X  :  float  :=  Low_X; 

Y  :  float  :=  Low_Y; 

Z  :  float  :=  Low_Z; 

Fofxyz  :  float; 

Pi  :  float  :=  3.141592654; 

begin 

for  k  in  1 . . (Count_Z)  loop 
Y  :=  Low_Y; 

for  j  in  1.. (Count _Y)  loop 
X  :=  Low_X; 

for  i  in  l..(Count_X)  loop 
if  (k  =  1)  then 

C_Netf(i,j,k)  :=  0.0; 

C_Calc(i, j ,k)  :=  0.0; 
elsif  (k  =  Count_Z+l)  then 
C_New(i, j ,k)  :=  0.0; 

C_Calc(i, j ,k)  :=  0.0; 
elsif  (j  *  1)  then 
C_New(i,j,k)  :*  0.0; 

C_Calc(i.j,k)  :=  0.0; 
elsif  (j  =  Count _Y+1)  then 
C_New(i, j ,k)  :=  0.0; 

C_Calc(i, j ,k)  :=  0.0; 
elsif  (i  =  1)  then 
C_New(i,j,k)  :=  0.0; 

C_Calc(i,j ,k)  :=  0.0; 
elsif  (i  =  Count_X+l)  then 
C_New(i, j ,k)  :=  0.0; 

C_Calc(i, j ,k)  :=  0.0; 
else 

if  (T  =  0.0  and  X  *  0.0  and  Y  <=  0.0 
and  (Y  +  Del.Y)  >0.0 

and  Z  <=  Height  and  (Z  +  Del_Z)  >  Height) 

then 

Fofxyz  :=  Quantity; 

else 

Fofxyz  :=  0.0; 
end  if; 
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Mu  :  float  :»  (A*Del_T)/(2.0*Del_X) ; 

Alpha  :  float  :=  (Kl*Del_T)/(Del_X*Del_X) ; 

Beta  :  float  :=  (K2*Del.T)/(Del_Y*Del_Y) ; 

Geunma  :  float  :=  (K3*Del_T)/(Del_Z*Del_Z) ; 

C_New(i,j,k)  :=  C_N(i,j,k)  -  Mu*(C_N(i+l, j ,k)  - 

C_N(i-l,j.k))  + 

Alpha*(C_N(i+l, j ,k)  -  2 .0*C^M(i , j ,k)  + 
C_N(i-l,j.k))  + 

Beta+(C_N(i, j+l,k)  -  2.0*C_N(i , j ,k)  + 
C_N(i,j-l,k))  + 

Gainnia*(C_N(i,j  ,k+l)  -  2.0*C_N(i,j  ,k)  + 
C_N(i,j  ,k-l))  + 

Fofxyz*Del_T ; 
if  T  /=  0.0  then 

C_Calc(i, j ,k)  :=  (Qu2Uitity/(8.0*sqrt(Kl*K2*K3) 
♦sqrt (Pi*Pi*Pi*T*T*T) ) ) * 
exp (- ( (X-A*T) ♦ (X-A*T) / (4 . 0*K1*T) )  - 
(Y*Y/(4.0*K2*T))  - 
(Z*Z/(4.0*K3*T))); 

end  if; 
end  if; 

X  :=  X  +  Del_X; 
end  loop; 

Y  :*  Y  +  Del.Y; 
end  loop; 

Z  :=  Z  +  Del_Z; 
end  loop; 

for  i  in  l..Count_X  loop 
for  j  in  l..Count_Y  loop 
for  k  in  1 . .Count_Z  loop 
C_N(i,j,k)  :=  C_New(i,j,k); 
end  loop; 
end  loop; 
end  loop; 

end  Compute_Final ; 


—  PROCEDURE:  Print.Info 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 

—  It  prints  out  the  matrix  with  the  points  for  the  time 

—  iteration  given. 

—  INPUT  PARAMETERS:  C_New  :  Tri_Matrix  for  calculating  answer 
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increment or  for  space 
incrementor  for  Y  direction 
incrementor  for  Z  direction 
incrementor  for  time 
lover  bound  for  X 
upper  bound  for  X 
lower  bound  for  Y 
upper  bound  for  Y 
lower  bound  for  Z 
upper  bound  for  Z 
Tri .Matrix  answer 
output  file  with  values  of  matrix 
Count.X,  Count.Y,  Count.Z  :  integer  counters 

X.Value,  Y.Value,  Z.Value  :  X,Y,and  Z  points 
Kl,  K2,  K3  :  diffusivity  coefficients 
BlMu,  B2Mu,  B3Mu  :  Stability  requirements  <*1/8 
Alphal,  Alpha2,  Alpha3  :  Qscilatory  condition  <*  1 
GLOBALS  USED:  None. 

CALLED  BY:  main. 

CALLS :  none . 


Del.X 

Del.Y 

Del.Z 

Del.T 

Low.X 

Hi_X 

Low.Y 

Hi_Y 

Low.Z 

Hi_Z 

OUTPUT  PARAMETERS:  C.New 

File 

LOCAL  VARIABLES: 


procedure  Print.Info (C.New  :  in  out  Tri.Matrix; 

C.Calc  :  in  out  Tri.Matrix; 

Del.X,  Del.Y,  Del.Z,  Del.T,  A  :  in  float; 
Low.X,  Hi.X,  Low.Y,  Hi.Y  :  in  float; 
Low.Z,  Hi.Z,  Time  :  in  float; 

Quantity,  Height  :  in  float; 

File  :  in  out  Text_Io.File.Type)  is 

Kl  :  float  :=  1.0; 

K2  :  float  :=  1.0; 

K3  :  float  :=  1.0; 

Count _X  :  integer  :=  integer((Hi_X  -  Low.X) /Del.X) ; 

Count _Y  :  integer  :=  integer ((Hi.Y  -  Low.Y) /Del.Y) ; 

Count.Z  :  integer  :=  integer((Hi_Z  -  Low.Z) /Del.Z) ; 

X. Value  :  float  :=  Low.X; 

Y. Value  :  float  :*  Low.Y; 

Z. Value  :  float  :=  Low.Z; 

BlMu  :  float  :=  (Kl*Del_T)/(Del_X*Del_X) ; 

B2Mu  :  float  :=  (K2*Del_T)/(Del_X*Del_X) ; 

B3Mu  :  float  :=  (K3*Del_T)/(Del_X*Del_X) ; 

Alphal  :  float  :=  (A*Del_X)/(2.0i'Kl) ; 
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Alpha2  :  float  :=  (A*Del_X)/(2 .0*K2) ; 
Alphas  :  float  :=  CA*Del_X)/(2.0*K3) ; 


begin 

Text.Io .put_line("#The  following  are  the  results  using  the 
forward  time  ") ; 

Text_Io.put_line("#  central  space  for  the  partial  diferential 
equation: ") ; 

Text_Io.put_line("#  Ct  +  A^Cx  =  Kl*Cxx  +  K2*Cyy  +  K3*C2Z  + 
F(X,Y,Z)"); 

Text_Io.put_line("#  The  resulting  matrix  is:  "); 

Text.Io . new.line ; 

Text_Io .new_line(File,2) ; 

Text.Io. put (File,"#  Output  for  threedprob.a  "); 

Text_Io .new_line(File) ; 

Text_Io.put_line(File,"#  Ct  +  A*Cx  =  Kl*Cxx  +  K2*Cyy  +  K3*Czz 
+  F(X.Y.Z)"); 

Text.Io .new_line(File) ; 

Text.Io. put (File,"#  Min_X  =  "); 
My_Float_Io,put(File,Low_X,3,2,0) ; 

Text.Io. put (File,",  Max_X  =  ") ; 

My_Float_Io . put (File ,Hi_X ,3,2,0); 

Text_Io. put (File,",  Delta.X  =  ") ; 
My_Float_Io.put(File,Del_X,3,5,0) ; 

Text_Io.new_line(File) ; 

Text. lo. put (File,"#  Min.Y  =  ") ; 

My.Float.Io .put(File,Low_Y,3,2,0) ; 

Text_Io.put(File,",  Max.Y  =  ") ; 

My.Float.Io .put (File,Hi_Y,3,2,0) ; 

Text _Io. put (File,",  Delta.Y  =  ") ; 

My.Float.Io .put (File, Del.Y, 3, 5,0) ; 

Text.Io .new.line(File) ; 

Text. lo. put (File,"#  Min.Z  =  ") ; 

My.Float.Io .put (File, Low.Z, 3, 2,0) ; 

Text.Io.  put  (File,",  MaJc.Z  =  ")  ; 

My.Float.Io . put (F ile ,Hi_Z ,3 ,2,0); 

Text.Io. put (File,",  Delta.Z  =  "); 

My.Float.Io. put (File, Del. Z, 3, 5,0) ; 

Text.Io .new.line (File) ; 

Text.Io. put (File,"#  Delta.T  =  "); 

My.Float.Io .put (File, Del.T, 3, 3,0) ; 

Text.lo.new.line(File) ; 
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Text_Io.put(File,"#  Time  of  calculation  =  ") ; 
My_Float_Io.put(File,Time,3,3,0) ; 

Text_Io.neH_line(File) ; 

Text. lo. put (File,"#  A  *  ") ; 

My_Float_Io.put(File,A,3,2,0) ; 

Text_Io.new_line(File) ; 

Text _Io. put (File,"#  Emission  Value  ®  "); 

My.Float.Io. put (File, Quantity, 3, 2,0) ; 
Text_Io.new_line(File) ; 

Text _Io. put (File,"#  Stack  height  =  "); 

My.Float.Io. put (File, Height ,3,2,0) ; 

Text. lo.new.line (File) ; 

Text. lo. put (File,"#  BlMu  =  (Kl*Del.T)/(Del_X*Del_X)  =  ") 
My.Float.Io.put(File,BlMu,3,4,0) ; 

Text.Io.new.line(File) ; 

Text. lo. put (File,"#  B2Mu  =  (K2*Del.T)/(Del.X*Del.X)  =  ") 
My. Float. lo. put (File, B2Mu, 3, 4,0) ; 

Text.Io.new.line(File) ; 

Text. lo. put (File,"#  B3Mu  =  (K3*Del.T)/(Del.X*Del.X)  =  ") 
My.Float.Io . put (File, B3Mu ,3,4,0); 

Text.Io.new.line(File) ; 

Text. lo. put (File,"#  Alphal  =  (A*Del.X)/(2*Kl)  »  ") ; 
My.Float.Io.put(File,Alphal,3,4,0) ; 

Text. lo.new.line (File) ; 

Text. lo. put (File,"#  Alpha2  =  (A*Del_X)/(2*K2)  =  "); 
My.Float.Io .put(File,Alpha2, 3, 4,0) ; 

Text. lo.new.line (File) ; 

Text. lo. put (File,"#  Alpha3  =  (A*Del_X)/(2*K3)  =  ") ; 
My.Float.Io .put(File,Alpha3, 3, 4,0) ; 

Text.Io.new.line(File) ; 

Text.Io.set.col(File,l) ; 

Text. lo. put (File,"#  X  "); 

Text.Io.set.col(File,10) ; 

Text. lo. put (File,"  Y  "); 

Text.Io.set.col(File,20) ; 

Text. lo. put (File,"  Z  "); 

Text.Io . set.col(File ,30) ; 

Text. lo .put (File, "Numerical") ; 

Text.Io. set .col(File, 47) ; 

Text.Io. put (File, "Actual") ; 

Text.Io.new.line(File) ; 
for  i  in  1.. Count _X+1  loop 
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Z.Value  :=  Low_Z; 
for  k  in  l..Count_Z+l  loop 
Y_ Value  :«  Low_Y; 
for  j  in  1.. Count _Y+1  loop 
Text_Io.8et_col(File,l) ; 

My_Float_Io. put (File, X_Value,l ,2,0) ; 
Text_Io.set_col(File,10) ; 

My_Float_Io .put(File,  Y_Value, 1 ,2 ,0) ; 
Text.Io . set_col(File,20) ; 

My_Float_Io. put (File,  Z_Value, 1,2,0) ; 
Text_Io.set_col(File,30) ; 

My_Float_Io. put (File, C_Netf(i,j ,k) ,1,6,2) ; 
Text_Io . set_col(File,47) ; 
My_Float_Io.put(File,C_Calc(i, j ,k) ,1,6,2) ; 
Text_Io.new_line(File) ; 

Y.Value  :*  Y.Value  +  Del.Y; 
end  loop; 

Z_Value  :=  Z_Valua  +  Del_Z; 
end  loop: 

X,Value  :=  X.Value  +  Del.X; 
end  loop; 

Text_Io.new_line; 
end  Print.Info; 


"  PROCEDURE:  Main 

—  DESCRIPTION:  This  is  the  main  progreun  emd  calls  all  the  above 

modules 

"  INPUT  PARAMETERS:  none 
"  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES:  none 
"  GLOBALS  USED:  All. 

—  CALLED  BY :  none 

—  CALLS:  Get.Info,  Compute_Final,  Print.Info. 


begin  — MAIN 

Text_Io . Create (Outfile ,Text_Io .0ut_File , "testSd .dat") ; 
Get_Info  (Delta.X,  Delta.Y,  Delta_Z,  Delta_T,  Min_X,  Max_X, 
Min.Y,  Max.Y,  Min_Z,  Max_Z,  A,  T.Value, 

Q_Value,  Stk_Hgt) ; 

T  :=  0.0; 

while  T  <=  T_Value  loop 

Compute_Final(C_N,  C_New,  C_Calc,  T,  Delta_X,  Delta_Y, 
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Delta_Z,  Delta_T,  Min_X,  M2Ut_X,  Min_Y,  Max_Y, 
Min.Z,  Mzut.Z,  A,  Q.Value,  Stk_Hgt) ; 

T  :=  T  +  Delta.!; 
end  loop; 

Print_Info(C_New,  C.Calc,  Delta.X,  Delta.Y,  Delta.Z,  Delta.!, 
A,  Min.X,  Max.X,Min.Y,Max.Y,Min_Z,Max.Z,!.Value, 
Q.Value,Stk_Hgt,Outfile) ; 

!ext.Io. Close (Outfile) ; 
end  threed; 
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B.5  Steady-State  Equation  Ada  Code 

—  FILE:  ubareqkcon.a 

—  PROJECT:  Thesis  work,  example  calculation 
"  DATE:  15  Sep  93 

—  VERSION:  Version  1.0 

—  AUTHORS :  Capt  Dave  Paal 

—  DESCRIPTION:  This  program  solves  a  partial  differential 

equation  using  the  forward  time  zmd  central  space. 

The  partial  differential  equation  to  be  solved  is: 

U_Bar*Cx  =  KY*Cyy  +  KZ*Czz 

where  KY/U_Bar  =  KZ/U_Bar  =  1.0 

and  (Ki*Del_T)/(U_Bar*Del_X*Del_X)<=  1/4 

"  with  C(0,x,y)  =  sin((Pi/20)*(Y  +  10))  *  sin((Pi/20)*(Z  +  10)) 
and  C(x,y,z)  =  0.0  at  Y=+-10  emd  Z=+-10 
"  OPERATING  SYSTEM:  UNIX/Sun  Sparc  Station 
"  LANGUAGE:  Meridian  Ada 
"  FILES  USED:  Output:  kconequ.dat. 


—  CONTEXT  CLAUSES 

with  text_io; 
with  Math.Lib; 
use  Math_Lib; 
with  my_integer_io; 
with  my_float_io; 

procedure  kconequbar  is 

—  TYPE  DECLARATION 

type  Matrix  is  array  (0. .90,0. .90)  of  float; 


—  GLOBAL  VARIABLES  AND  EXCEPTIONS 
C_N  :  Matrix; 

C_New  :  Matrix; 

C_Calc  :  Matrix; 

X  :  integer; 
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Delta_X 

float ; 

Delta_Y 

float ; 

Delta_Z 

float ; 

Min.X 

float ; 

Max_X 

float ; 

Min.Y 

float ; 

Max_Y 

float ; 

Min_Z 

float ; 

Max_Z 

float ; 

U.Bar 

float ; 

Outf ile 

Text_Io .File_Type; 

—  PROCEDURE:  Get  information  for  calculation 

—  DESCRIPTION:  This  procedure  gets  the  values  of  Del_X,  Del_Y, 

Del_Z,  Low_X,  Hi_X,  Low_Y,  Hi_Y,  Low_Z,  Hi_Z,  U_Bar, 
Quauitity,  Height. 

—  INPUT  PARAMETERS:  None. 


— 

OUTPUT  PARAMETERS:  Del.X 

:  incrementor  for  space 

— 

Del.Y 

:  incrementor  for  Y  direction 

— 

Del.Z 

:  incrementor  for  Z  direction 

— 

Low_X 

;  lower  bound  for  X 

— 

Hi_X 

:  upper  bound  for  X 

— 

Low.Y 

:  lower  bound  for  Y 

— 

Hi_Y 

:  upper  bound  for  Y 

— 

Low.Z 

:  lower  bound  for  Z 

— 

Hi.Z 

:  upper  bound  for  Z 

— 

Ubar 

:  constant  wind  speed 

— 

LOCAL  VARIABLES :  None . 

—  GLOBALS  USED:  None. 

—  CALLED  BY :  main . 

—  CALLS :  none . 


procedure  Get_Info  (Del_X  :  in  out  float; 
Del_Y  :  in  out  float; 

Del_Z  :  in  out  float; 

Low_X  :  in  out  float; 

Hi_X  :  in  out  float; 

Low_Y  :  in  out  float; 

Hi_Y  :  in  out  float; 

Low_Z  :  in  out  float; 

Hi_Z  :  in  out  float; 

Ubar  :  in  out  float)  is 
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begin 


Text_Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

X.  "); 

Hy_Float_Io . get (Del_X) ; 

Text.Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

Y.  •'); 

My_Float_Io.get (Del_Y) ; 

Text_Io .put_line("Enter  the  value  (Floating  point)  of  Delta 

Z.  "); 

My_Float_Io.get(Del_Z) ; 

Text_Io .put_line("Enter  the  minimum  value  (Floating  point) 
of  X.  "); 

My_Float_Io.get (Low_X) ; 

Text_Io .put_line("Enter  the  value  (Floating  point)  of  X 
wanted.  ") ; 

My_Float_Io.get (Hi_X) ; 

Text_Io.put_line("Enter  the  minimum  value  (Floating  point) 
of  Y.  "); 

My_Float_Io.get(Low_Y) ; 

Text_Io.put_line("Enter  the  maximum  value  (Floating  point) 
of  Y.  ") ; 

My_Float_Io.get(Hi_Y) ; 

Text.Io .put_line("Enter  the  minimum  value  (Floating  point) 
of  Z.  "); 

My_Float_Io . get (Low_Z) ; 

Text _ I o. put .line ("Enter  the  maximum  value  (Floating  point) 
of  Z.  "); 

My_Float_Io .get (Hi_Z) ; 

Text.Io .put_line("Enter  the  value  (Floating  point)  of  U_Bar. 

"); 

My_Float_Io.get(Ubar) ; 
end  Get.Info; 

—  PROCEDURE:  Compute  final  Matrix 

—  DESCRIPTION:  This  procedure  calculates  final  answer  through 

iterations  using  a  finite  difference  scheme  with  forward 

—  time/central  space. 

—  INPUT  PARAMETERS:  C_N  :  Matrix  for  C_N(i,j) 

C_New  :  Matrix  for  calculating  numerical 
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C_Calc  :  Matrix  for  calculating  exact 
answer 

Del_X  :  incrementor  for  space 
Del_Y  :  incrementor  for  Y  direction 
Del_Z  :  incrementor  for  Z  direction 
Low_X  :  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 

where  calculation  is  made 
Low_Y  :  lower  bound  for  Y 
Hi_Y  :  upper  bound  for  Y 
Low_Z  :  lower  bound  for  Z 
Hi_Z  :  upper  bound  for  Z 
Ubar  :  constant  wind  speed 

—  OUTPUT  PARAMETERS;  C_N  :  Matrix  for  C(N) 

C_New  :  Matrix  for  calculating  numerical 
amswer 

C_Calr.  :  Matrix  for  calculating  exact 
answer 

—  LOCAL  VARIABLES:  KY  :  Ubar*sigmaY*2/2*X 

KZ  :  Ubar*sigmaZ“2/2*X 

Alpha  :  (KY*Del.X)/(Ubar*Del_Y*Del_Y) 
Beta  :  (KZ*Del.X)/(Ubar*Del.Z*Del_Z) 
Count.Y  :  integer((Hi_Y  “Low_Y)/Del_Y) 
Count_X  :  integer ((Hi_X  -  Low_X)/Del_X) 
Count.Z  :  integer ((Hi_Z  -  Low_Z)/Del_Z) 

X  :  value  of  x 

Y  :  value  of  y 

Z  :  value  of  z 

Pi  :  the  number  Pi 

—  GLCBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 


procedure  Compute_Final(C_N  :  in  out  Matrix; 

C_New  :  in  out  Matrix; 

C_Calc  :  in  out  Matrix; 

X  :  in  integer; 

Del.X,  Del_Y,  Del.Z  :  in  float; 
Low_X,  Hi_X,  Low_Y,  Hi_Y  :  in  float; 
Low_Z,  Hi_Z  :  in  float; 

Ubar  :  in  float)  is 


KY  :  float  :=  Ubar; 

KZ  :  float  :=  Ubar; 

Count.X  ;  integer  :*  integer((Hi_X  -  Low_X)/Del_X) ; 

Count_Y  :  integer  :=  integer((Hi_Y  -  Low_Y)/Del_Y) ; 

Count_Z  ;  integer  :®  integer((Hi_Z  -  Low_Z)/Del_Z) ; 

X_Value  :  float  :=  Low_X  +  float (X-l)*Del_X; 

Y  :  float  :=  Low_Y; 

Z  ;  float  :*  Low_Z; 

Alpha  :  float  :=  (KY*Del_X)/(Ubar*Del_Y*Del_Y) ; 

Beta  :  float  :=  (KZ*Del_X)/(Ubar*Del_Z*Del_Z) ; 

Fofxyz  :  float; 

Pi  :  float  :=  3.141592654; 

PiDiv20  :  float  :=  Pi/20.0; 

begin 

Z  :=  Low_Z; 

for  k  in  1 . . (Count_Z)  loop 
Y  :=  Low.Y; 

for  j  in  1 . . (Count.Y)  loop 
if  (k  =  1)  then 

C_New(j ,k)  :=  0.0; 
elsif  (k  =  Count_Z+l)  then 
C_New(j ,k)  :=  0.0; 
elsif  (j  =  1)  then 
C_New(j ,k)  :=  0.0; 
elsif  (j  =  Count_Y+l)  then 
C_New(j,k)  :=  0.0; 
elsif  (X.Value  =  0.0)  then 
C_New(j,k)  := 

sin((PiDiv20)*(Y+10.0)''*sin((PiDiv20)*(Z+10.0))  ; 
else 

Alpha  :  float  :=  (KY*Del_X)/(Ubar*Del_Y*Del_Y) ; 

Beta  :  float  :=  (KZ*Del_X)/(Ubar*Del_Z*Del_Z) ; 
C_New(j,k)  :=  C_N(j,k)  + 

Alpha*(C_N(j+l,k)  -  2.0*C.N(j,k)  + 
C_N(j-l,k))  + 

Beta*(C_N(j ,k+l)  -  2.0*C_N(j,k)  + 
C_N(j,k-l)); 

end  if; 

C_Calc(j,k)  :=  exp((-(Pi*Pi)*X_Value)/20))* 

sin( (PiDiv20) ♦ (Y+10 . 0) ) ♦sinC (PiDiv20) ♦ (Z+10 . 0) ) ; 
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Y  :=  Y  +  Del.Y; 
end  loop; 

Z  Z  +  Del_Z; 
end  loop; 

for  k  in  l..(Count_Z)  loop 
for  j  in  1 . . (Count_Y)  loop 
C_N  C j , k)  ; =  C_New ( j , k) ; 
end  loop; 
end  loop; 

end  Compute_Final; 


—  PROCEDURE:  Print.Info 

—  DESCRIPTION:  This  procedure  prints  the  output  of  the  program. 

—  It  prints  out  the  matrix  with  the  points  for  the  time 

—  iteration  given. 

—  INF”  PARAMETERS:  C_New  :  Matrix  for  calculating  numerical 

einswer 

C_Calc  ;  Matrix  for  calculating  exact 

—  answer 

Del.X  :  incrementor  for  space 
Del_Y  :  incrementor  for  Y  direction 
Del.Z  :  incrementor  for  Z  direction 

—  U_Bar  :  advection  constemt 

Low_X  ;  lower  bound  for  X 
Hi_X  :  upper  bound  for  X 
Low_Y  :  lower  bound  for  Y 
Hi_Y  :  upper  bound  for  Y 
Low_Z  :  lower  bound  for  Z 
Hi_Z  :  upper  bound  for  Z 

—  OUTPUT  PARAMETERS:  C_New  :  Matrix  numerical  aoiswer 

C_Calc  :  Matrix  for  calculating  exact 
emswer 

File  :  output  file  has  values  of  matrix 

—  LOCAL  VARIABLES:  Count_X,  Count _Y,  Count _Z  :  integer  counters 

X.Value,  Y_Value,  Z.Value  :  X,Y,  and  Z  points 

—  GLOBALS  USED:  None. 

—  CALLED  BY:  main. 

—  CALLS :  none . 

procedure  Print.Info (C_New  :  in  out  Matrix; 

C.Calc  :  in  out  Matrix; 
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Del.X,  Del_Y,  Del_Z,  Ubar  :  in  float; 

Low_X,Hi_X,Low_Y,Hi_Y,Low_Z,Hi_Z  :  in  float; 
File  :  in  out  Text_Io.File_Type)  is 
KY  :  float  :*  Ubar; 

KZ  :  float  :=  Ubar; 

Count_X  :  integer  :=  integer((Hi_X  -  Lou_X)/Del_X) ; 

Count _Y  :  integer  :=  integer((Hi_Y  -  Lou_Y)/Del_Y) ; 

Count.Z  :  intege*"  :=  integer((Hi_Z  -  Low_Z)/Del_Z) ; 

X. Value  :  float  :=  Low_X; 

Y. Value  :  float  :=  Low_Y; 

Z. Value  :  float  :=  Low_Z; 

StableY  :  float  :=  (KY*Del_X)/(Ubar*Del_Y*Del_Y) ; 

StableZ  :  float  :=  (KZ*Del.X)/(Ubar*Del_Z*Del_Z) ; 


begin 

Text_Io.put_line("#The  following  are  the  results  using  the 
forward  time  ") ; 

Text_Io.put_line("#  central  space  for  the  partial  diferential 
equation: ") ; 

Text_Io.put_line("#  Ubar*Cx  =  KY*Cyy  +  KZ*Czz  ") ; 

Text_Io.put_line("#  The  resulting  matrix  is:  ") ; 
Text_Io.new_line; 

Text_Io.new_line(File.,2) ; 

Text  _Io.  put  (File,''#  Output  from  file  ubareqkcon.a  ")  ; 
Text_Io.new_line(File) ; 

Text_Io.put_line(File,"#  Ubar*Cx  =  KY*Cyy  +  KZ*Czz  ") ; 
Text_Io.new_line(File) ; 

Text_Io.put_line(File, "#  KY  2uid  KZ  are  constant  emd  equal 
U_bar") ; 

Text.Io .new_line(File) ; 

Text.Io. put (File,"#  Min_X  =  "); 

My_Float _Io . put (File , Low_X ,3,2,0); 

Text _Io. put (File,",  Max_X  =  ") ; 

My_Float_Io .put(File,Hi_X,3,2,0) ; 

Text.Io. put (File,",  Delta.X  =  "); 

My_Float_Io .put (File,Del_X,3,5,0) ; 

Text_Io.new_line(File) ; 

Text.Io. put (File,"#  Min.Y  =  ") ; 

My .Float. lo .put(File,Low.Y,3,2,0) ; 

Text.Io. put (File,",  Max.Y  =  ") ; 

My.Float.Io .put(File,Hi.Y,3,2,0) ; 

Text.Io. put (File,",  Delta.Y  =  ") ; 
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My_Float_Io .put(File,Del_Y,3,5,0) ; 

Text_Io .new_line(File) ; 

Text.Io. put (File,"#  Min_Z="); 
My_Float_Io.put(File,Low_Z,3,2,0) ; 

Text _Io. put (File,",  Max_Z  =  ") ; 
My_Float_Io.put(File,Hi_Z,3,2,0) ; 
Text_Io.put(File,",  Delta_Z  *  ") ; 

My_Float_Io .put (File,Del_Z,3,5,0) ; 
Text_Io.new_line(File) ; 

Text_Io. put (File,"#  U_Bar  =  "); 

My_Float_Io . put (File ,Ubar ,3 , 5,0); 
Text_Io.new_line(File) ; 

Text_Io .put (File, "#  Y  stability  value  =  "); 
My_Float_Io.put(File,StableY,3,5,0) ; 
Text_Io.new_line(File) ; 

Text_Io .put (File, "#  Z  stability  value  ®  "); 
My_Float_Io .put(File,StableZ,3,5,0) ; 
Text_Io.new_line(File) ; 

Text_Io.set_col(File,l) ; 

Text.Io. put (File,"#  X  "); 

Text_Io.set_col(File,10) ; 

Text.Io. put (File,"  Y  "); 

Text.Io. set.col (File, 20) ; 

Text.Io. put (File,"  Z  ") ; 

Text.Io.set.col(File,30) ; 

Text.Io. put (File, "Numerical") ; 

Text.Io . set_col(File,47) ; 

Text.Io. put (File, "Actual") ; 

Text.Io .new_line(File) ; 

for  k  in  1.. Count .Z+1  loop 
Y.Value  :=  Low.Y; 
for  j  in  1.. Count .Y+1  loop 
Text.Io. set .col(File,l) ; 

My .Float. lo. put (File, Hi.X,l ,2,0) ; 
Text.Io. set .col(File, 10) ; 

My.Float.Io .put (File,  Y.Value, 1,2,0) ; 
Text.Io . set. col(File, 20) ; 

My.Float.Io .put(File,  Z.Value, 1 ,2 ,0) ; 
Text.Io. set. col(File, 30) ; 

My.Float.Io. put(File,C.New(j ,k) ,1,6,2) ; 
Text.Io .  set.coKFile, 47)  ; 

My.Float.Io .put (File, C.Calc(j ,k) , 1 ,6,2) ; 
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Text _Io .new_line (File) ; 
Y.Value  :=  Y_ Value  +  Del.Y; 
end  loop; 

Z_Value  :*  Z_Value  +  Del_Z; 
end  loop; 

Text_Io.new_line; 
end  Print_Info; 


—  PROCEDURE:  Examplel .3. 1  (a.k.a  main) 

—  DESCRIPTION:  This  is  the  main  progrzua  and  calls  all  the  above 

modules 

—  INPUT  PARAMETERS:  none 

—  OUTPUT  PARAMETERS:  none 

—  LOCAL  VARIABLES:  none 

—  GLOBALS  USED:  All. 

—  CALLED  BY:  none 

—  CALLS:  Get_Info,  Compute.Final,  Print.Info. 


begin  — MAIN 

Text_Io.Create(Outfile,Text_Io.Out_File, "kconequ.dat") ; 
Get_Info  (Delta.X,  Delta.Y,  Delta.Z,  Min.X,  Mm.X,  Min^Y, 
Max.Y,  Min.Z,  Max.Z,  U.Bar) ; 

X  :=  integer ((Max_X  -  Min_X)/Delta_X)+l; 
for  I  in  1 . .X  loop 

Compute_Final(C_N,  C_New,  C.Calc,  I,  Delta.X,  Delta.Y, 

Delta_Z,  Min_X,  M2uc_X,  Min.Y,  Max_Y,  Min_Z,  Hax_Z,  U_Bar) ; 
end  loop; 

Print_Info(C_New,  C.Calc,  Delta.X,  Delta_Y,  Delta.Z,  U.Bar, 

Min.X,  Max_X,  Min_Y,  Max.Y,  Min.Z,  Max.Z,  Outf ile) ; 
Text.Io.Close(Outfile) ; 
end  kconequbar; 
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Appendix  C.  1-D  diffusion  and  3-D  steady  state  equations 

This  appendix  describes  the  one-dimensional  diffusion  equation  and  the  three- 
dimensional  steady  state  equations  and  how  they  used  to  this  research. 

C.l  The  Diffusion  Equation: 

The  one-dimensional  diffusion  equation  (C.l),  also  known  as  the  one-dimen¬ 
sional  heat  equation,  is  the  simplest  parabolic  partial  differential  equation. 


Ct  =  bcx 


(C.l) 


The  finite  difference  scheme  used  in  this  research  to  solve  equation  (C.l)  is  the 
forward-time  central-space  scheme 


«n+l  _  -n 

At 


r”  —  9/’”  4-  r" 

I  1 

Ax2 


(C.2) 


which  can  also  be  written  as 


c:+‘ = c + -  2c + c_. )  (C.3) 

where  //  =  AtjAx^.  The  stability  condition  for  this  method  is  kxfi  <  1/2.  It  is 
found  by  replacing  with  in  Equation  C.2  leaving 


C-I 


g-l  e*^  -  2  +  e-'^ 
At  ~  ^  Ax2 

and  solving  for  g  gives 


(C.4) 


g  =  I  —  Akxfi  sin^  ^6.  (C-5) 

The  condition  <  1  from  Theorem  2.2.1  in  Strikwerda  (26:42)  is  equivalent  to 
4:kxfi  sin^  <  2  which  is  true  for  all  6  if  and  only  if  kifi  <  1  /2.  Thus  this  method 
is  conditionally  stable. 

C.2  3-D  Steady  State  Equation 

The  Eulerian  approach  (25:542)  to  the  three-  dimensional  steady  state  equation 
examined  in  this  research  is 

UCx  -  kxCxX  4"  kyCyy  4"  kgCyy  “j-  ^  (C.6) 

with 

c(x,t,,z)  =  0  x,y,z^±oo  (C.7) 

where  u  is  the  wind  speed  (positive  constant),  kx  —  ky  =  k,,  are  constant  in  this 
research  and  are  the  x-axis,  y-axis,  and  z-axis  dillusivity  terms  respectively,  and  q  is 
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the  source  term.  The  forward-time  central-space  finite  difference  scheme  (C.8)  used 
to  solve  equation  (C.6)  is 


The  solution  to  eq.  (C.8)  would  be  compared  to  the  exact  solution  (25:542)  of 


q  .  u(r-x), 

c  =  __ —  expl - — -] 

Airkir  2ki 


(C.9) 


where  q  —  the  source  term,  ki  —  (kxkykz)^^^,  jind  r  =  {x^  + 


The  stability  condition  is  the  similar  to  the  three-dimensional  advection-diffu- 
sion  discussed  in  Section  3.4.4.  It  is  assumed  that  k^Cxx  C  ucx,  that  is,  the  diffusion 
in  the  x-axis  direction  is  much  less  than  the  advection  in  the  x-  axis  direction  and 
is  therefore  negligible.  It  is  for  this  reason  that  this  equation  was  not  looked  at  in 
detail. 
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